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A  normally  incident  Rayleigh  wave  may  be  used  for  the  investigation  of  a  generauL 
vertical  boundary  and  the  attenuation  of  a  viscoelastic  medium.  Use  of  energy 
conservation  emd  proper  boundary  conditions  produce  2  N  second  order 
differential  equations,  N  being  the  number  of  viscoelastic  layers  in  the  medium. 
The  homogeneous  part  of  the  differential  equations  can  be  transformed  into  an 
eigenvalue  problem  by  the  use  of  finite  element  technique;  the  eigenvalue  amd 
eigenvectors  of  the  eigenvalue  problem  are  the  wavenumbers  and  the  displacement 
amplitudes  of  the  viscoelastic  la 


red  medium.  The  real  and  imaginary  parts 
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the  wavenumber  determine  the  phase  velocity  and  the  attenuation  of 
the  layered  medium,  respectively.  Dispersion  and  the  attenuation  curves 
can  be  obtained  by  using  different  periods.  The  above  wavenumbers  can  be 
used  in  the  inhomogeneous  differential  equation;  this  equation  contains 
the  effect  of  the  vertical  boundary;  solution  of  this  equation  determines 
the  displacement  at  the  vertical  boundary.^  The  displacements  at  the  left 
and  right  sides  of  the  boundary  may  be  used  in  the  energy  equation  to 
determine  the  reflected  and  the  transmitted  energy,  respectively.  The 
difference  between  the  incident  and  transmitted  energy  determines  the 
difference  between  the  wavenumbers  at  the  left  and  right  hamd  sides  of 
the  vertical  boundary.  Addition  of  these  differences  to  the 
wavenumbers  of  the  left  hand  boundary  would  result  in  the  formation  of 
•  the  dispersion  curve  at  the  right  hand  boundary. 
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ABSTRACT 


'  A  normally  incident  Rayleigh  wave  may  be  used  for  the 
investigation  of  a  general  vertical  boundary  and  the  attenua- 
tion  of  a  viscoelastic  medium.  Use  of  energy  conservation  and 
proper  boundary  conditions  produce  2  N  second  order  differen¬ 
tial  equations,  N  being  the  number  of  viscoelastic  layers  in 
the  medium.  The  homogeneous  part  of  the  differential  equations 
can  be  transformed  into  an  eigenvalue  problem  by  the  use  of 
finite  element  technique;  the  eigenvalue  and  eigenvectors  of 
the  eigenvalue  problem  are  the  wavenumbers  and  the  displacement 
amplitudes  of  the  viscoelastic  layered  medium.  The  real  and 
imaginary  parts  of  the  wavenumber  determine  the  phase  velocity 
and  the  attenuation  of  the  layered  medium,  respectively.  Dis¬ 
persion  and  the  attenuation  curves  can  be  obtained  by  using 
different  periods.  The  above  wavenumbers  can  be  used  in  the 
inhomogeneous  differential  equation;  this  equation  contains 
the  effect  of  the  vertical  boundary;  solution  of  this  equation 
determines  the  displacement  at  the  vertical  boundary.  The 
displacements  at  the  left  and  right  sides  of  the  boundary  may 
be  used  in  the  energy  equation  to  determine  the  reflected  and 
the  transmitted  energy,  respectively.  The  difference  between 
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the  incident  and  transmitted  energy  determines  the  difference 
between  the  wavenumbers  at  the  left  and  right  hand  sides  of 
the  vertical  boundary.  Addition  of  these  differences  to  the 
wavenumbers  of  the  left  hand  boundary  would  result  in  the 
formation  of  the  dispersion  curve  at  the  right  hand  boundary. 
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PREFACE 


The  attenuation  of  a  surface  wave  along  its  path  is  of  interest  in 
many  branches  of  physical  science;  among  the  noted  are  solid  mechanics, 
acoustics,  and  seismology.  The  aim  of  this  dissertation  is  to  present 
a  method  of  Rayleigh  wave  propagation  in  a  viscoelastic  solid  con¬ 
strained  by  one  or  two  vertical  boundaries.  The  material  is  arranged 
to  present  a  suninary  of  viscoelastic  behavior  of  solid,  an  introduction 
to  the  method  of  finite  element  for  the  derivation  and  solution  of 
governing  equations,  a  discussion  on  energy  conservation  in  elastic  and 
viscoelastic  solids,  and  presentation  of  a  perturbation  theory  for  the 
calculation  of  the  phase  velocity  due  to  a  vertical  boundary,  followed 
by  a  chapter  on  application  of  data  to  the  method. 

I  hope  that  this  dissertation  will  be  a  stimulant  for  interested 
investigators,  to  be  critical  of  the  material  herein,  and  to  be  encour¬ 
aged  enough  to  further  test  the  validity  of  the  theory  by  using  other 
distinct  models.  As  a  reference,  the  text  is  supolemented  by  a  rela¬ 
tively  vast  number  of  references,  which  provide  a  great  deal  of 
information  for  understanding  the  subject. 

I  would  like  to  thank  Professor  E.  Herrin  for  his  advice  through 
the  course  of  my  research  and  for  his  critical  opinion  of  this  dis¬ 
sertation.  I  would  also  like  to  thank  Doctors  B.  Mohraz,  D.  Blackwell, 
T.  Goforth,  H.  Mack,  and  W,  Peeples,  for  their  criticisms  and 
encouragements. 
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I  also  wish  to  express  my  appreciation  to  Dr,  L.  Drake;  his 
computer  program  and  his  personal  assistance  greatly  contributed  to  my 
research.  Finally,  I  am  grateful  to  my  wife,  Jeneva,  for  her  patience 
and  encouraging  attitude  during  this  study. 

Support  for  this  study  was  provided  by  the  Defense  Advanced 
Research  Projects  Agency  under  Contract  AFOSR  F49620-76-C-0030  monitored 
by  the  U.S.  Air  Force  Office  of  Scientific  Research;  by  National  Science 
Foundation  Grant  EAR-77-13714;  and  by  Institute  for  the  Study  of  Earth 
and  Man  Grant  20-91-235. 
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CHAPTER  I 


INTRODUCTION 

Rayleigh  wave  amplitude  decays  as  it  propagates  through  the  earth. 
The  amount  of  the  decay  observed  is  invariably  more  than  what  the  theory 
predicts.  This  anomaly  is  due  to  the  loss  mechanisms  or  the  energy 
dissipation  in  the  medium;  thus  the  observed  Rayleigh  wave  amplitude 
is  a  function  of  energy  dissipation  or  attenuation. 

Attenuation  of  Rayleigh  wave  is  due  to  anelasticity  and  geo¬ 
metrical  discontinuities  of  the  earth.  Anelasticity  of  the  earth  may 
be  thought  of  as  the  imperfection  in  crystal  structure  such  as  point 
defects,  dislocation  across  which  the  atoms  are  not  aligned  in  accor¬ 
dance  with  the  normal  lattice  structure,  and  interatomic  bonds  at  grain 
boundaries.  Geometrical  discontinuities,  among  others  (23,  24),  cause 
the  distortion  of  the  plane  of  propagation  of  Rayleigh  waves  from  the 
great  circle  and  therefore  there  is  a  loss  of  energy  or  attenuation. 

In  a  dissipating  medium  the  deformation  of  a  solid  is  a  function 
of  time,  temperature  and  space.  Creep  phenomenon  is  a  deformation  process 
in  which  there  is  a  time  dependent  strain  under  a  constant  applied  stress 
and  relaxation  phenomenon  is  a  time  dependent  stress  under  a  constant 
strain.  A  solid  manifesting  creep  and  relaxation  phenomena  is  called 
viscoelastic.  If  a  deformation  is  due  to  discontinuities  in  the  medium 
or  the  amplitude  decay  is  only  a  function  of  position,  the  dissipation 
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of  the  wave  energy  1s  of  the  hysteretic  type  (Figure  1),  and  the  amount 
of  this  loss  must  be  determined  by  the  average  energy  dissipation  In  a 
hysteresis  loop  over  one  cycle  of  harmonic  motion  and  Is  a  function  of 
displacement  amplitude. 

To  narrow  the  gap  between  the  observed  and  predicted  Rayleigh 
wave  amplitude  the  effect  of  viscoelasticity  and  geometrical  discon¬ 
tinuities  must  be  taken  Into  account.  A  few  workers  have  studies  the 
viscoelastic  behavior  of  Rayleigh  waves  for  an  Isotropic  homogeneous 
layered  half-space.  Press  and  Mealy  (94)  derived  an  expression  for  the 
absorption  coefficient  and  a  dispersion  relation  assuming  complex 
velocities  for  shear,  longitudinal  and  Rayleigh  waves  respectively. 

King  and  Sheard  (56)  used  a  dissipation  function  for  a  general  anlso- 
trop'.c  medium  and  derived  a  dispersion  relation  for  an  Isotropic  layered 
half-space  equivalent  to  that  of  Press  and  Mealy  (94).  Borcherdt  (15) 
also  formulated  a  similar  theoretical  expression  for  a  layered  half¬ 
space.  Auld  (8)  determined  the  attenuation  of  Rayleigh  waves  due  to 
viscous  damping  In  the  medium  by  using  King  and  Shread’s  dissipation 
function  (56).  He  concluded  that  the  Rayleigh  wave  attenuation  Is 
mostly  due  to  shear  viscosity  of  the  medium.  For  observational  studies 
of  the  dissipation  of  Rayleigh  waves  the  reader  Is  referred  to  Anderson, 
et  al  (5)  and  Tryggvason  (115). 

The  effect  of  geometrical  discontinuities  on  the  propagation  of 
Rayleigh  waves  have  been  a  problem  of  seismologists  for  many  years. 
Several  studies  have  been  devoted  to  the  description  of  the  field  around 
wedges  and  corners  using  approximate  analytical  techniques:  e.g.,  Lapwood 
(66),  Mai  and  Knopoff  (82),  McGarr  (85)  and  Alsop  et  al  (2).  Several 
theoretical  studies  have  also  been  done  In  more  recent  years 


with,  in  most  part,  the  assumption  of  weak  horizontal  Inhomogeneity  on 
assumption  of  a  small  change  In  elastic  parameters.  Among 
the  more  Important  are.  Babish  et  a1.  (9),  Kennett  (55)  and  Hukitina 
et  al.  (88). 

The  observed  distorted  Rayleigh  wave  through  the  continental 
margins  or  the  oceanic  ridges  has  been  studied  by  Capon  (19,  20)  and 
by  Evemden  (31  ,  32,  33)  and  the  amount  of  distortion  has  been  approx¬ 
imated  analytically  by  Gjevik  (39)  who  has  used  ray  tracing  theory  and 
perturbation  theory.  The  distortion  of  Rayleigh  waves  from  the  great 
circle  has  also  been  studied  by  Crampin  (23,  24)  who  has  stated  that 
the  distortion  Is  due  to  a  layer  of  anisotropic  olivine  In  the  upper 
mantle. 

The  finite  element  method  has  been  used.  In  recent  years.  In 
wave  propagation  problems  of  seismology.  Lysmer  (78)  and  Lysmer  and 
Waas  (79)  developed  a  lumped  mass  method  to  study  the  propagation  of 
Rayleigh  and  shear  waves  1n  a  homogeneous  layered  medium.  Lysmer  and 
Drake  (80)  introduced  the  method  of  consistent  mass  matrix  of  the  finite 
element  method  to  seismology.  Drake  (28,  29)  studied  the  propagation 
of  Love  and  Rayleigh  waves  in  a  nonhorirontally  layered  medium.  Drake 
(26)  also  studied  the  motion  of  Rayleigh  waves  at  a  continental 
boundary.  Waas  (117)  Improved  the  method  of  finite  element  to  a 
general  two  dimensional  problem  and  Kausel  et  al.  (54)  extended  the 
method  to  a  three  dimensional  problem  In  which  Rayleigh  wave  corres¬ 
ponds  to  a  vertical  cross  section  and  Love  wave  corresponds  to  a 
horizontal  cross  section  of  the  medium.  Smith  (107)  applied  the  method 
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of  finite  element  to  body  wave  propagation.  Finally  Segol  et  al. 
(105)  studied  the  propagation  of  a  Rayleigh  wave  being  obstructed  by 
trenches . 

In  this  dissertation  the  finite  element  method  is  used  to  study 
the  attenuation  of  Rayleigh  waves  due  to  viscoelasticity  i.e.,  vis¬ 
cosity,  temperature,  partial  melting,  etc.,  and  determine  a  dissipative 
factor  due  to  a  vertical  boundary.  The  model  consists  of  a  heterogeneous 
medium  bounded  by  two  vertical  boundaries  (Figure  4).  These  two  bound¬ 
aries  are  the  common  boundaries  between  the  heterogeneous  medium  and 
two  semi  infinite  horizontal  layered  media  extending  to  left  and  right. 

For  the  two  viscoelastic  layered  media  a  second  order  differential 
equation  similar  to  that  of  a  damping  oscillator  is  obtained  in  which 
the  damping  factor  is  a  matrix  whose  entries  are  the  viscoelastic  para¬ 
meters  of  the  medium.  The  solution  of  this  equation  determines  the  wave- 
numbers  and  the  amount  of  attenuation  of  Rayleigh  waves  in  the  horizontal 
direction.  Use  of  the  boundary  conditions  at  the  vertical  boundaries 
results  in  the  equation  of  motion  of  the  heterogeneous  medium  whose 
solution  yields  the  horizontal  and  vertical  displacement  of  the 
heterogeneous  medium.  These  displacements  and  the  tvavenumbers  of  the 
layered  regions  are  used  in  the  energy  conservation  theorem  along  with 
a  perturbation  theory  to  determine  the  amount  of  absorption  due  to 
viscoelasticity  of  the  region  and  the  change  in  the  wavenumbers  due  to 
vertical  boundaries.  The  adjusted  wavenumbers  of  the  right  hand 
boundary  may  be  used  to  determine  the  dispersion  curve  and  the  amount 
of  attenuation  at  the  right  hand  boundary. 


CHAPTER  II 


VISCOELASTICITY 
A.  Suimiary 

The  phenomenon  of  viscoelasticity  Is  a  broad  subject  In  the  study 
of  materials  In  the  gaseous,  liquid,  and  solid  states.  Interested 
readers  are  referred  to  the  selected  texts  on  the  subject  by  Bland  (13), 
Christensen  (22),  Gross  (44),  Kolsky  (60,  61),  Litvitz  and  Davis  (72), 
Mathson  (84)  and  Zener  (122)  for  theoretical  studies,  and  Anderson  (3, 
4),  Anderson  et  al.  (5,  6),  Gordon  (41,  42,  43),  Hart  et  al.  (45), 
Jackson  and  Anderson  (50),  Liu  et  al.  (74),  Magnitsky  and  Zharkov  (81) 
and  Orowan  (90)  for  seismological  Interest.  A  brief  review  for  the 
problem  at  hand  will  be  presented  here. 

B.  Uniaxial  Deformation 

Viscoelastic  materials  combine  elastic  and  viscous  effects; 
therefore.  It  is  possible  to  model  the  material  of  the  viscoelastic 
earth  by  networks  made  of  simple  elastic  and  viscous  elements.  The 
simpllest  model  Is  a  spring  and  a  dashpot  In  series  or  parallel  and 
Is  called  a  Maxwell  or  Voigt  element,  respectively.  Mathematically 
the  stress-strain  relation  for  the  model  Is 


for  the  Maxwell  element,  and 
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we  +  n^c  »  o  (2.2) 

for  the  Voigt  element,  where  o  is  the  stress,  e  is  the  strain,  c  the  rate  of 
the  strain,  w  is  the  elastic  modulus  of  the  spring  and  n  is  the  coefficient 
of  viscosity  of  the  dashpot  (Figure  2).  A  combination  of  a  spring  with 
elastic  modulus  w^  either  in  series  with  a  Voigt  element  or  in  parallel 
with  a  Maxwell  element  will  produce  a  relation 

+  o  *  M^(c  +  t^c)  (2.3) 

called  the  standard  linear  solid,  where 


w  *  ”r  w^'^  w 


(2.4) 


are  the  stress  relaxation  time  under  constant  strain,  the  strain 
relaxation  under  constant  stress,  and  the  deformation  or  relaxed  elastic 
modulus,  respectively. 

The  solution  of  Equation  2.1  for  a  constant  strain  results  in  a 
decaying  exponential  function  of  the  form 

a(t)  -  Ojj  EXP  (-  w/nyt)  (2.5) 


which  is  called  the  relaxation  function.  Similarly  the  solution  of  Equation 
2.2  for  a  constant  stress  results  in  an  increasing  exponential  function 
of  the  form 

”  w/n^t 

e(t)  »  (1  -  e  M  (2.6) 

which  is  called  a  creep  function  (Figure  3).  The  solution  of  Equation  2.3 
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yields  the  relaxation  and  creep  functions  of  the  form 

a(t)  -  M_c(t)  +  M  ^EXP  (-  t/T  )  f\x?  (t/t  )E(T)dT 

»  I  T  O  O 

O  •• 

and  (2.7) 

E(t)  »  jjp  o(t)  +  ^  ^  EXP  (-  t/tp)  /  EXP(t/t  )o(T)dT 

Equations  2.7  represent  the  Boltzmann  after-effect.  The  verbal  statement 
of  these  equations  by  Boltzmann  Is  that  for  all  real  solids  the  relation 
between  stress  or  strain  components  at  any  time  depends  not  only  on  their 
Instantaneous  values  at  that  time*  but  on  the  whole  history  of  defor¬ 
mation  from  the  time  at  which  the  material  was  formed.  It  is  now 
possible  to  represent  the  viscoelastic  behavior  of  the  real  material 
by  either  a  large  number  of  Maxwell  elements  and  springs  joined  In 
parallel  or  a  large  number  of  Voigt  elements  and  springs  joined  In  series. 
The  arrangement  of  the  springs  represent  the  molecular  position  of 
material  and  the  dashpots  represent  the  perturbation  factor  for  the 
distortion  of  molecular  positions.  Thus  the  generalized  relaxation  and 
creep  function  for  uniaxial  deformation  Is  given  by 

N  t 

a(t)  -  M^c(t)  +  Mj.  r  EXP  (-  t/r^  )  /  EXP(t/t^  )c(T)dT 
1*1  1  “•  1 

and  (2.8) 

N  ^  t 

c(t)  «  a(t)  +  tJ-  r  ~  exp  (-  t/T^  )  /  EXP  (t/t^  )o(T)dT 
"r  ”r  1«1  S'  «  S 

Equations  2.8  represent  the  Boltzmann  superposition  principle. 


C.  Three  Dimensional  Deformation 


For  a  general  linear  viscoelastic  model  consider  the  energy  density 

(67), 

*  “  ^ijKL  *^10. 

where  c's  are  strains  and  generally  complex,  are  the  moduli  of  the 

medium.  For  an  anisotropic  medium  contain  21  independent  complex 
coefficients.  For  an  isotropic  medium  these  coefficients  reduce  to  two 
elastic  constants  in  an  elastic  medium,  whereas,  in  the  viscoelastic 
medium  the  two  coefficients  are  not  constant  and  may  be  a  function  of 
time,  temperature,  position  or  other  thermodynamic  parameters.  There¬ 
fore,  in  an  isotropic  viscoelastic  medium  Equation  2.9  can  be  written  as 

V/  ■  1/2  A  ^i1^ 

+  w'(t,  T,  i,  (2.10) 

where  A  and  v  are  Lame*  constants.  A*  and  v'  the  coefficients  which 
depend  on  the  viscoelastic  behavior  of  the  medium.  Since  Equation  2.10 
is  a  scalar  function,  the  stress  field  can  be  obtained  by  a  directional 
derivative  with  respect  to  the  strain  tensor  as 

(T^j  *  +  2  uc^j  +  A  ( t,  T,  X,  +  2  u  (t,T,x,Y)c^j 

(2.11) 

The  first  two  terms  in  Equation  2.11  represent  the  familiar  form  of  the 
stress-strain  relation  in  an  isotropic  homogeneous  medium  and  the  last 
two  terms  are  due  to  the  viscoelasticity  of  the  medium. 
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D.  Effective  Moduli  for  a  Relaxing  Medium 

For  functionals  x'  and  w*,  the  effective  moduli  of  the  stressed 
medium  can  be  obtained  by  either  substitution  of  Equation  2.11  into 
Equation  2.7  and  integrating  it  through  the  limits  of  integration,  or 
by  formation  of  a  relation  similar  to  Equation  2.3,  using  Equation  2.11; 
the  latter  is  a  more  convenient  way  to  choose.  For  this  purpose  we  may 
write  Equation  2.11  in  terms  of  shear  and  dilatational  components. 

Using  the  strain  tensor 

"ij  “  ("ij  •  ^  .  (2.12) 

Equation  2.11  becomes 

°ij  *  ^ij^ii^  *  *ij^ii^ 

+  (2.13) 

where  K  ■  x  +  2/3ii  is  called  dilational  or  bulk  modulus.  Separating 
the  stress  tensor  of  Equation  2.13  into  the  compressional  and  sheer 
component,  the  result,  in  the  form  of  Equation  2.3,  can  be  written  as 

•  • 

°ii  *  %  "ii  ‘  Hi  -  %  "ii 

(2.14) 

®ij  ^  "s  °ij  “  ‘'"ij  *  S  "ij 

where  t  and  t.  are  relaxation  times,  K  «  K  and  v  *  u^v/t.  .  The 

assumption  of  a  harmonic  stress  in  Equations  2.14  will  produce  a 
relation  of  the  form 
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o^^(a))  *  -  Kg((i))  c^^(u) 

(2.15) 

a^j(a))  »  Ug(w)  c^j(w) 

(2.16) 

where 

Kg(«)  -  K[1  -  (1  - 

-  1u  (1  -  ^S/t^)  '^v/(1  +  W^T^^)] 

pg(<*i)  *  p  [1  -  (1  -  )i<)^Tj^/(1  +  w^Tg^) 

-  iw  (1  -  %/Tj)Tj/(1  +  W^Tg^)] 

(2.17) 

Equation  2.17  is  identical  to  Equation  2.16  of  Liu  et  al.  (74) 

which  was  derived  Independently.  It  is  clear  now  that  the  effective 

moduli  of  a  viscoelastic  wave  and  therefore  the  velocity  of  the  wave  is 

a  function  of  frequency  which  has  been  claimed  by  many  investigators 

for  body  waves  and  recently  by  a  few  for  surface  waves.  Futterroan  (38) 

and  Lomintz  (77)  theorized  the  concept  of  frequency  dependence  of  the 

phase  velocity  of  body  waves.  Kogan  (59),  Lairt>  (65),  Liu  et  al.  (74), 

Magnitsky  and  Zharkov  (81),  Mason  (83)  and  Strick  (109)  have  studied  the 
frequency  dependent  velocity  of  body  waves  and  concluded  the  validity  of 
Futterman's  theory.  Most  recently  a  few  studies  have  been  devoted  to 
frequency  dependent  phase  velocity  of  surface  wave  and  the  free  oscil¬ 
lation  of  the  earth  by  Liu  and  Archambeau  (73.  75),  Park  and  Rockwell 
(92),  Randall  (95)  and  Kanamori  and  Anderson  (53). 
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E.  Two  Dimensional  Deformation 


The  six  independent  components  of  the  stress  tensor  in  a  three 
dimensional,  isotropic  medium  reduce  to  three  independent  components 
in  two  dimensions.  Therefore,  an  x  -  z  plane  propagation  would  result 
in  a  zero  stress  in  y  direction  such  that 


xy  yy 


zy 


(2.18) 


Starting  with  the  first  two  terms  of  Equation  2.11  (for  derivation  of 
the  equations  of  motion,  the  inelastic  part  of  the  stress  tensor  need 
not  be  included  for  the  moment). 


°ij  '  "ii  ^ 


i  =  1,3 
j  *  1,3 


(2.19) 


and  the  constraint  relation.  Equation  2.18,  the  nonzero  three  indepen¬ 
dent  component  of  the  stress  in  two  dimensional  systems  are 


’ll 


'XX 


+  2w  +  Xe 


XX 


XX 


'rz 


^  *  o  *  2u  c 
13  xz  xz 


‘’33  *  ‘’zz  '  ^"xx  ^  2*'  ‘=zz"^^'=zz 


(2.20) 


where  the  components  of  the  strain  tensor  are  obtained  from  the  relation 

(2.21) 


3U.  3U. 


i  »  1,3 
j  *  1,3 


as 


*4K>S*V 


3U^ 

l/2(^ 

+  — 2-) 
3X  ’ 

3U, 

3U, 

3U. 

dU 

‘■11  ^xx 


■13  ■  ^X2 


^33  =  ^22  '  -  3i^)  =  35^ 


(2.22) 


where 

Ux  and  are  hori2ontal  and  vertical  displacements  of  the  propagating 
plane  wave  respectively.  Substituting  Equations  2.22  into  Equations  2.20 
and  expressing  the  stress  components  in  a  vector  form,  the  result  is 


°  =  (^X*  °22‘  °XZ^  =  ^  " 


(2.23) 


where 


X+2u  X  0 

.  /3U^  3U,  3U,  3U  \ 

0=  X  X.2u  0 

_  0  0  p  _ 


(2.24) 


Equation  2.23  is  derived  for  an  elastic  medium;  for  an  inelastic  medium 
the  Lame'  constants  in  Equation  2.23  must  be  replaced  by  the  effective 
moduli  of  Equation  2.17. 

F.  Wave  Equation  in  Cartesian  Coordinates 

Using  the  stress  tensor  of  Equation  2.19  the  wave  equation  is 
obtained  by 


’U.J  ■  ‘“u’  ‘ 


(2.25) 


If  the  force  F^.  is  assumed  to  be  the  inertia  force  then 


T -  ( O  ■  •  )  -  pU . 

3Xj  '  ir  ^  1 


(2.26) 


where  p  is  the  density  and  double  dot  denotes  the  acceleration. 
Assumption  of  a  hormonic  plane  wave  allows  a  displacement  of  the  form 


where 


U.  =  U.(x,z)e^“^ 


U.  =  (U^.U^) 


(2.27) 


Substitution  of  Equations  2.20,  2.22  and  2.27  into  Equation  2.26  leads  to 
a  two  dimensional  equation  of  motion 


2. 1 

-  ou  u. 


(2.28) 


Since  the  displacement  is  a  function  of  x  and  z,  an  independent 
product  of  the  function  of  x  and  z  may  be  replaced  by  U^. (x,z) 

Ux(x,z)  =  U^(z)W(x) 

(2.29) 

U2(’'»z)  =  U2(2)W(x) 

Substitution  of  Equations  2.29  into  Equations  2.28  yields  the  coupled 
ordinary  second  order  differential  equations 
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u 


?  2 
iK  (X  +  m)  ^  -  r(X  +  2u)Uj^  + 


0 


(x  +  2u) 


2  2 

iK  (X  +  u)  ^  -  KSU^  +  p/u^  =  0 


(2.30) 


and 


^  +  K^W 
dx^ 


0 


Equation  2.31  has  a  solution  of  the  form 
W(x)  = 


(2.31) 


(2.32) 


and  Equation  2.27  may  be  expressed  by 
=  U^(2)e^<“^  ■ 

=  U^(z)e^^“^  ■ 


(2.33) 


the  amplitudes  U  (z)  and  U  (z)  are  obtained  by  the  solution  of  the 
^  £ 

differential  Equations  2.30. 

So  far  there  has  not  been  any  restriction  on  the  condition  of 
Equation  2.30  and  the  amplitudes  U^(-)  and  U^(z)  Equation  2.30  may  be 
solved  for  a  general  plane  wave.  Our  interest  is  In  a  solution  with 
proper  constraints  to  produce  a  Rayleigh  wave,  a  plane  wave  that  propa¬ 
gates  In  one  direction  and  dissipates  In  the  direction  normal  to  the 
direction  of  propagation.  In  general,  these  two  directions  need  not  be 
perpendicular;  when  the  medium  of  propagation  Is  viscoelastic  the 
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direction  of  propagation  and  dissipation  of  Rayleigh  wave  are  not 
perpendicular.  In  this  case  the  wave  Is  called  a  generalized  Rayleigh 
wave  and  Is  discussed  In  more  detail  In  Chapter  IV. 

To  acquire  a  Rayleigh  wave  solution.  In  an  N  layered  medium,  from 
Equation  2.33,  we  analyze  Equation  2.30.  For  an  N  layered  medium  Equation 
2.30  constitutes  2N  second  order  differential  equations.  The  boundary 
conditions  required  for  the  solution  of  these  equations  are  continuity 
of  the  stresses  and  displacements  at  each  Interfaces,  zero  stress  at 
free  surface,  and  zero  displacement  amplitude  at  a  distance  far 
from  the  surface.  The  stresses  considered  are  the  stress  normal  to  the 
plane  Interface  and  the  shear  stress  on  the  plane  interface;  they  are 
given  by 


0..  »  X  ^  +  (X  +  2u)  ^ 


zz 


3Z 


xz 


(2.34) 


There  must  be  continuity  at  the  Interfaces  where  the  Lame'  constants  just 
above  and  below  the  Interfaces  may  be  different.  Thus,  the  conditions  for 
the  stresses  are 


w* 


an  f  ^ 


on  the  surface  and 

3U 


3U.  3U^ 

^j-1  3x^  ^^j-1  ^“j-l^  iT  “  ^’j-1  W  *  ^^'j-1  ^  ’'j-1'  3Z 


-1  \  3Z  ^  3X  z  ^j-1  '  3Z  3X  / 


3U  3U. 


(2.36) 
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j  =  2,...  N,  at  the  interfaces.  Also  the  conditions  for  continuity 
of  displacements  are 

U,  «  U'  and  U,  *  Ul  j  «  2....  N  (2.37) 

j  j  j  j 

and  for  far  enough  below  the  surface 

U,  -  0  and  U,  «  0  .  (2.38) 

*N+1  ^N+1 

With  these  boundary  conditions  the  solution  of  Equation  2.30  will  lead  to 
an  eigenvalue  problem  whose  solutions  are  extremely  difficult.  For  this 
reason  we  leave  further  discussion  of  Equation  2.30  for  the  next  chapter 
in  which  we  discuss  the  finite  element  method. 


I 


1 


CHAPTER  III 


THEOREHCAL  BACKGROUND  OF  THE  FINITE  ELEMENT  METHOD 
A.  Sunnary 

For  a  detailed  review  of  the  method  of  finite  element  the  reader 
Is  referred  to  the  texts  by  Bathe  and  Wilson  (11),  Desal  and  Abel  (25), 
Norrie  and  de  Vries  (89),  Richards  (97),  Segerlind  (103)  and  Zlenklewlcz 
(124)  for  mathematical  and  engineering  application,  and  by  Lysmer  (78), 
Lysmer  and  Waas  (79),  Lysmer  and  Drake  (80),  Segol  et  a1,  (105)  and 
Waas  (117)  for  selsmologlcal  applications.  A  sunmary  of  the  derivation 
of  the  essential  equations  will  be  discussed  in  this  chapter. 

The  finite  element  technique  has  been  brought  to  seismology  in 
recent  years.  It  seems  to  be  more  successful  than  the  method  of  finite 
difference.  The  basic  difference  between  the  two  techniques  is  that 
in  the  finite  difference  method  the  continuous  medium  Is  used  to  derive 
a  set  of  differential  equations  which  are  replaced  by  finite  difference 
equations.  In  the  finite  element  technique  the  continuous  medium  is 
first  approximated  and  the  governing  differential  equations  are  solved 
numerically.  Therefore,  In  the  finite  element  technique  the  basic 
assumption  Is  to  divide  the  continuous  layered  medium  into  a  discrete 
system  consisting  of,  for  example,  rectangular  elements  each  having 
four  nodal  points.  The  displacement  of  the  structure  is  equivalent  to 
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the  displacement  of  the  nodal  points  and  all  the  forces  between  the 
elements  and  the  external  forces  are  transmitted  through  the  nodal 
points. 

The  forces  that  produce  displacements  within  an  element  are 
represented  by 

F  -  KU  (3.1) 

where  l(  Is  called  the  stiffness  matrix  depending  on  the  orientation  and 
properties  of  the  elements.  If  the  Inertia  force  Is  present  Equation  3.1 
Is  written  as 

F  «  (K  -  M)  U  (3.2) 

where  M  Is  the  mass  matrix  of  the  element  and  u  Is  the  frequency  of  the 
vibration.  Since  the  number  of  elements  have  to  be  finite  the  structure 
must  be  bounded.  For  this  reason  the  structure  with  Irregular  geometry, 
I,  Is  joined,  along  the  vertical  planes,  to  the  regions,  L  and  R,  which 
extend  Infinitely  to  the  left  and  to  the  right,  respectively  (Figure  4). 
Regions  L  and  R  consist  of  horizontal  layers  which  are  welded  together 
at  their  Interfaces  and  may  differ  In  their  material  properties  and 
thickness.  For  this  configuration  the  finite  element  method  Is  used 
to  analyze  the  motion  In  the  Irregular  zone  I  and  a  discrete  theory, 
which  Is  essentially  a  one  dimensional  finite  element.  Is  used  to 
analyze  the  motion  of  the  layered  zones.  The  natural  layers  of  the  L 
and  R  zones  must  be  subdivided  Into  a  larger  number  of  layers  so  that 
the  nodal  points  of  the  elements  in  the  Irregular  zone  at  the  vertical 
boundaries  coincide  with  the  layer  Interfaces.  The  thickness  of  the 
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element  Is  determined  by  the  wave  length  of  the  S  wave  in  the  layered 
zones  (104). 


B.  Setup  of  the  Method 

The  Irregular  region  of  zone  I  consists  of  an  assemblage  of 
quadrilateral  elements  interconnected  at  their  nodal  points.  A  typical 
element  e  is  defined  by  its  nodes  i,  j,  K  and  L.  The  displacements  at 
the  element  nodes  are  taken  as  unknown  and  the  solution  within  each 
element  is  approximated  by  a  shape  function  N  that  varies  linearly  along 
the  boundaries,  i.e.. 


G  =  N V  *  [N.,  Nj,  N^]  (V..  V^,  V^)  (3.3) 

in  which,  for  a  plane  wave  propagation, 

-  (V^.  V2)i  etc.  (3.4) 

are  the  horizontal  and  vertical  displacements  at  the  element's  node  i, 
etc.,  and  the  raartix  is  the  shape  function  in  terms  of  (e  •  n),  the 
local  coordinates  of  the  element.  The  component  N^,  etc.  of  the  matrix 
is  unity  at  the  node  i,  etc.  and  is  equal  to  zero  at  other  nodes. 

The  above  displacement  is  related  to  the  strain  by 


e 


BU 


where 


U  »  (u^.u^) 


(3.5) 


(3.6) 


\ 
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,If  the  strain  energy,  or  the  energy  density  ,  defined 
in  Chapter  I IB,  exists,  then  as  in  Chapter  I IB,  its  derivative  with 
respect  to  the  strain  tensor  is  the  stress  field  and  is  given  by 


3Cij  “ij 


(3.7) 


The  elastic  part  of  the  stress  in  Equation  3.7  in  terms  of  the  strain 
is  given  by 


°ii  '  ^"ij  ^  "ij 
Equation  3.7  can  be  written  equivalently  as 


(3.8) 


w  -  (3.9) 

which  is  the  virtual  work  done  by  the  force  per  unit  of  volume  and 
the  surface  force  per  unit  of  area  (37).  For  static  equilibrium  the 

expression 


/  6  c^jdv  *  /  F^  «  U^.dv  +  /  6  U^ds  (3.10) 

is  called  the  principle  of  virtual  work  produced  by  the  virtual  dis¬ 
placement  6U^.  Equation  3.10  will  be  used  to  derive  the  equation  of 
motion  of  an  element  in  the  irregular  region  I.  If  the  body  force  F^  is 
the  inertia  force  per  unit  of  volume, a  set  of  inhomogeneous  differential 
equations  can  be  derived  by  rewriting  Equation  3.10  as 


^  ^ij  "  ^i  ®  (3.11) 

This  equation,  with  T^  equal  to  zero  can  also  be  used  to  derive  a  homo¬ 
geneous  second  order  differential  equation  for  the  layered  medium,  L  and 
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R  of  Figure  4,  as 


/(o..  4  e^.  -  F.  U.)dv 


«  0  . 


(3.12) 


The  first  term  of  this  equation  can  be  written  as 


/  6  e^jdv  =  !  c..(6U.^  *  6Uj.)dv 


®  T  iU.dv  +  /  o. .  V  61). ds, 

^  I J  »J  1  I J  I 


(3.13) 


using  the  familiar  expression 


j  3X..  / 


J  1 


(3.14) 


and  the  Gauss ' s  theorem.  Substitution  of  Equation  3.13  into  Equation  2.12 


gives 


/(Cij^j  +  F.)  6U^dv  -  /  V  6U^dS  =  0 


(3.15) 


which  for  an  arbitrary  6U^  satisfies  the  equation  of  motion 


»ij.j  *  ^  ^ 


(3.16) 


and  either 


6U^  »  0 


(3.17) 


on  the  rigid  boundary  or 


"(j  '■  ■ 


(3.18) 


the  surface  boundary  condition. 
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Equation  3.15  determines  the  motion  of  the  wave  in  the  horizontal 
layers  L  and  R.  The  first  integrand  of  this  equation  is  the  general 
equation  of  motion  and  the  second  integrand  is  the  horizontal  boundary 
condition,  i.e.,  in  a  plane  wave  case,  the  normal  stress  in  the  z 
direction  and  the  shear  stress  on  the  plane  interfaces.  For  convenience 
we  write  Equation  3.15  in  a  vector  form  as 

/  .  W  dv  -  /  6U^  •  o  C  ds  =  0  (3.19) 

v  s 

where  W  is  the  equation  of  motion,  and  T  represents  complex  conjugate 
transpose. 


C.  Rayleigh  Wave  Motion 


The  displacement  of  the  medium  in  the  horizontal  and  vertical 
directions  produced  by  the  propagation  of  a  Rayleigh  wave  can  be 
represented  by  (1) 


“x  '  ‘’x'"-^! 

(3.20) 

Substitution  of  this  equation  into  the  equations  of  motion,  with  no  force 
present,  and  the  separation  of  variables  In  the  usual  way  would  result  in 


Uj^(x,z)  «  Uj^(z)  e’^*^ 

U2(x,z)  •  U2(z)  e'^*^ 


(3.21) 


Substitution  of  Equation  3.21  into  Equation  3.20  will  produce  the  vector 
displacement  of  Equation  3.19  as 
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U  •  (Ujj(2).  U^( 


z))e 


-  Kx) 


(3.22) 


The  six  component  stress  tensor  in  three  dimensions  will  reduce  to  a 
three  component  stress  vector  in  two  dimensions  and  is  given  by 


°  *  <^x*  “zz*  °xz) 


where,  from  Equation  3.8 


XX 


3U  3U, 
«  (X  +  2u)  r-^  +  X 


3X 


3U 


az 

au. 


°zz  “  ^  af  "  ^ 


/au  au  V 


a  * 
XZ 


Therefore, 


0*0 


where 


x+2tj 

X 

0 


X  0 

X+2y  0 
0  u 


.('!j<  !lli 

'  ax  ’  3z  ’ax  3z  ' 


(3.23) 


(3.24) 


(3.25) 


(3.25a) 


(3.26) 


At  the  interfaces  the  stress  component  in  the  x  direction  is  zero  and 
the  surface  stress  av  in  Equation  3.19,  and  the  matrix  D  in  this  case  are 


I 
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ov 


(o 


zz'  °xz 


) . 


(3.27) 


"  X  x+2p 
-0  0 


then 


ov  *  D' E 


(3.28) 


(3.29) 


respectively.  Substituting  Equations  3.22,  3.24-3.26  and  3.27-  3.29  Into 
Equation  3.19  and  assuming  that  Is  the  Inertia  force  per  unit  of 
volume  we  have 

ol  -  fij  58^  •  U)dxd2  -  /  D*  e  dx  (3.30) 

V  s 


where  the  elemental  volume  dv  and  surface  ds  In  three  dimension 
reduce  to  dxdz  and  dx  In  two  dimensions,  respectively.  Using  Gauss's 
theorem  on  the  last  term  of  Equation  3.30  we  get 

/  6U^*  O'  e  dx  «  /  DZ^  («U^«  D*  e  )dxdz  (3.31) 

s  v 

T  3  P  oT  3 

where  DZ  “  ^  J  represents  transpose.  Since  operates 

on  a  product  vectors  <U  and  c  the  right  hand  side  of  Equation  3.31 
expands  to 

/  Iy  ^(SU^- J'e  )dxdz  «  /  DZ^  (D'e  )dxdz 

V  " 


+  /  O'  DZ^  r  dxdz  (3.32) 

V 
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If  we  rewrite  Equation  3.6 


3U,  aU  \ 

'ax  az  / 


0 


az 

ax 


(u 


(3.33) 


-T 

c  * 


(U 


X’ 


(3.34) 


and  substitute  Equations  3.32  -  3.34  into  Equation  3.30  we  have 


/(«U^-  H  U  -  DZ’’’  U^*H'  U)  dxdz  *  0 


where 


H  = 


K^(a  +  2u)  ■  pii)^ 


3l 


„2  2 

K  U  -  Pu 


and 


H' 


iKX 

az 


(A  2w)  - 


iKp 


(3.35) 


(3.36) 


(3.37) 


Equations  3.35  involves  partial  derivatives  with  respect  to  z  of  the 
displacement  1)  in  Equation  3.22  which  is  an  implicit  function  of  z.  The 
differential  equations  obtained  in  this  way  will  possess  the  unknowns 


26 


in  the  argument  of  the  exponential  functions  and  thus  their  solution 
are  nearly  impossible.  For  this  reason  we  use  the  technique  of  the 
finite  element  to  derive  equations  that  may  be  solved  conveniently. 

For  the  layered  region  we  use  a  one  dimensional  element,  i.e., 
an  element  with  two  nodal  points  for  the  two  horizontal  plane  interfaces 
sharing  the  layer  as  shown  in  Figure  5.  Let  the  displacement 
vector  U  be  approximated  by, 

U  »  N  V  (3.38) 

so  that 

D’’’  -  (3.39) 

The  partial  derivative  of  this  with  respect  to  z  is 


3z 


where,  from  Equation  3.3,  for  a  linear  expansion 


N  =  [N.,  Nj]  =  1/2 


0  1  +  n 

-  n  0 


0 


(3.40) 


(3.41) 


If  the  origin  is  at  the  midpoint  of  the  layer  then 


h/2 

2  »  -h/2,  n  »  -1  at  the  top  of  the  layer  and  z  =  h/2,  n  »  1  at  the  bottom  of 
the  layer  (124)  as  shown  in  Figure  5. 
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1 


Substitution  of  Equations  3.38-3.40  into  Equation  3.35  results  in 


/fiV^(/KVa  N  +  iKN^  ~  N  -  pwVn  +  |^  H 


Layers 


+  iK  |r  N^V)  V  dxdz  =  0 

32  “ 


(3.43) 


where 


■x  +  2ii  0" 

"0  x" 

‘u  0  ■ 

i.  ® 

*  S-  * 

,  and  a  = 

_  0  y  _ 

o_ 

0  X  +  2y 

L  -J 

(3.44) 


Evaluation  of  the  integral  for  each  layer  gives 


/It  !l\  H  N •  £3 


h  T 

/  N  a  N  dz  *  Aj 


(3.45) 


"  T 

j  n'  N  dz  *  M . 
0 

And  from  Equation  3.42 


3_  ,  1 _ 3_ 

3z  h/2  3n 


(3.46) 


so  that 


—  N  *  ^  —  1/2 

3Z^  h/2  3n 


1  “n  0  1  +ri  0 


r  , 

-1  0  1  0 

V 1 

/ 

t] 

.0  -1  0  1. 

(3.47) 


1 


f 

'2(X+2w) 

0 

X+2v 

o“ 

'o 

X 

0  -x“ 

h. 

f  A  * 

:  "j  6 

0 

2u 

0 

V 

,  C.  =  1/2 

V 

0 

0 

A+2p 

0 

2(X+2m) 

0 

0 

X 

0  -X 

0 

p 

0 

2uj 

_  u 

0 

-1)  0_ 

Similarly, 


(3.49) 


=  p .  h  ./6 


The  layer  matrice  A.,  C.  and  M.  are  the  same  matrices  introduced 

U  V  V  w 

by  Lysmer  and  Drake  (80)  which  were  derived  originally  by  Lysmer  and 
Waas  (79)  by  a  different  method.  Addition  of  these  submatrices  over  the 


2  0  10 
0  2  0  1 
10  2  0 
0  10  2 


(3.50) 
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layers  by  the  method  of  Lysmer  and  Drake  (80),  as  shown  In  Figure  6 
for  a  typical  matrix  A.,  will  result  in 

J 

/  6V^(K^A  -  ilX  +  iKC^  +  6  -  «^M)  V  dx  «  0  (3.51) 

X 

For  an  arbitrary  virtual  displacement  vector  6V^  (37)  we  may  write, 
since  the  expression  in  parentheses  is  independent  of  x, 

(K^A  -  iKC  +  iKC^  +  G  -  A)  V  =  0  (3.52) 


Equation  3.52  represents  a  nonlinear  eigenvalue  problem  with  eigenvalue 
2 

K  and  eigenvector  V  whose  solution  is  usually  obtained  by  the  vector 
iteration  method.  The  analysis  of  the  derivation  of  the  solution  of 
Equation  3.52  will  not  be  presented  here  but  the  reader  is  recommended 
to  refer  to  the  excellent  texts  by  Bathe  and  Wilson  (11),  Waas  (117),  and 
Zienkiewicz  (124)  for  engineering  applications  and  by  Ruhe  (100), 
Papaconstantinov  (91),  and  Wilkinson  (120)  for  theoretical  discussions. 
At  a  particular  frequency,  the  generalized  eigenvalue  problem  (Equation 

3.52)  has  4N  eigenvalues  l  and  corresponding  eigenvectors  i  V., 

N  being  the  number  of  layers.  The  eigenvalues  may  be,  real,  imaginary 
or  complex. 

In  the  derivation  of  our  eigenvalue  problem  we  assumed  that  the 
matrix  moduli  0  be  real.  For  a  viscoelastic  model  we  may  substitute 
the  effective  moduli  of  Equation  2.17  of  Chapter  II  into  matrix  D, 
which  is 


X  +2m 

e  e 


^e*2we 


0 

0 


(3.53) 


0 


0 


where 

“e 

and  subscript  e  stands  for  effective  or  complex  moduli.  If  matrix  D  is 
real,  the  motion  of  the  Rayleigh  wave  in  the  layers  is  determined  by 
the  number  of  real  wave  numbers  obtained  from  the  solution  of  the  eigen¬ 
value  problem.  If  matrix  D  is  complex  the  eigenvalue  problem  has  complex 
eigenvalues.  The  real  part  of  the  eigenvalue  determines  the  phase 
velocity  of  the  wave  in  the  layers  and  the  imaginary  part  determines  the 
attenuation  of  the  wave  in  the  horizontal  direction.  Thus, 

phase  velocity  »  REAL°"twavenumber) 


attenuation  «  2 


I HAG  (wavenumber) 
real  ( wa  venumbe  r ) 


(3.55a} 


According  to  the  theory  of  linear  differential  equation,  if  each 
eigenvalue  and  the  corresponding  independent  eigenvectors  determine 
a  solution  then  a  linear  combination  of  these  solutions  is  also  a  solution. 
Therefore,  if  K.  for  1  *  1,  2N  are  the  eigenvalues  and  for  1  »  1,  2N 
are  the  corresponding  eigenvectors  the  complete  motion  of  the  layered 
region  may  be  obtained  by  superposition  of  these  solutions  as 


2N 
z  c 
i 


i  M 


gi(<ot  -  K^x) 


(3.56) 


in  which  each  multiplier  a  is  an  amplitude  factor  specifying  the 
contribution  of  the  mode  i  to  the  total  motion.  At  the  boundary  with 
no  phase  change  the  motion  is 
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where  V  Is  the  matrix  whose  columns  are  the  eigenvectors  V^.  If  the 
displacement  U  in  Equation  3.57  is  known. the  multiplier  vector  a  can 
be  solved  for  as 


o  =  V“^  U 


0.  Vertical  Boundary 


Let  us  rewrite  Equation  3.11 


/(o,j  Je,j  -  F,  6U,)dv  -  /T,  «U,d5 


and  substitute  Equation  3.13  into  it;  the  result  is 


/(Oij^j  +  F.)  «U^.dv  +  /(T^.  -  6U^.ds  =  0 

This  equation  can  be  satisfied,  for  an  arbitrary  if 


°fj.J  ♦  ■  0  . 

the  equation  of  motion,  and  either 


fiU.  «  0 


on  the  rigid  boundary  or 


/(T^  -  ^^i  *  0  » 


(3.58) 


(3.59) 


(3.60) 


(3.61) 


(3.62) 


(3.63) 


the  surface  boundary  condition  on  the  vertical  plane.  The  elemental 
surface  ds  in  this  case  is  dz,  therefore. 


/T^  6U^d2  =  fiU^dz 


(3.64) 


On  the  vertical  plane  the  stress  component  in  the  z  direction  is  zero 
and  the  surface  stress  o. .v  in  Equation  3.64,  in  vector  form,  is 

I J 

*  ^°xx*  °xz^ 

and  the  matrix  moduli  corresponding  to  Equation  3.28,  in  this  case  is 
r  X+2  X  0"1 

0^  *  (3.66) 

—  L  0  0  uJ 

Now  the  right  hand  side  of  Equation  3.63  becomes 

Jov  •  5U  dz  ■  D  e  dz  (3.67) 

s  s  — 

Substituting  Equations  3.33,  3.38-3.42,  and  3.57  into  Equation  3.67,  we 
have 

/ov  .  6U  dz  *  (i  A  V  K  -  C  V)  a  6V  (3.68) 

Layers  * 

The  derivation  of  this  equation  takes  the  same  steps  that  lead  to  the 
derivation  of  Equation  3.52.  In  Equation  3.68  all  matrices  and  the 
vector  a  have  been  defined  except  for  the  matrix  K.  Matrix  K  is  a 
diagonal  matrix  whose  diagonal  elements  are  the  eigenvalues  obtained 
from  the  solution  of  the  eigenvalue  problem.  Equation  3.68  represents 
the  force  applied  at  the  boundary  with  the  assumption  that  the  finite 
element  region  is  removed.  The  reaction  forces  acting  on  the  finite 
element  region  follow  then  from  the  equilibrium  consideration,  i.e., 
a  force  equal  and  opposite  to  that  of  Equation  3.68 
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Substitution  of  a  from  Equation  3.58  Into  Equation  3.69  yields 

P  =  L  U  (3.70) 

where 

L  -  1  A  V  K  -  C  (3.71) 

Is  called  the  stiffness  matrix  of  the  vertical  boundary. 


E.  Motion  In  the  Finite  Element  Region 

Rewriting  Equation  3.10  with  the  assumption  that  the  Inertia 
force  be  ,  we  have.  In  vector  form, 

/  o  •  6c  dv  -  U)  dxdz  =  /6U^.  Tv  dz  (3.72) 

V  s 


The  left  hand  side  of  Equation  3.72  can  be  written  as 


/(6c^*  De  -  PO)  6U^«  U)  dxdz  «  /(6uV«  0  B  U  -  pa.^  U)  dxdz 
V  V  “ 

(3.73) 

where  matrix  £  Is  given  by  Equation  3.6.  Substituting  Equations 
3.38-3.42  Into  the  expression  3.73,  Equation  3.72  becomes 


/(6yT^  0  B  N^V  -  pj-  6v\  •  N^V)  dxdz  =  +  /6v\  •  Tv  dz 


(3.74) 


where,  from  Equation  3.73,  for  a  two  dimensional  expansion 


~*sn(n-l)  0  l-n^  0  >sn(ri+l)  0  "1 

N  «  CN.,N.,N„,N,3  .  2 

^  ^  L  0  *5n{n-l)  0  l-n^  0  %n(n+l)J 


(3.75) 
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The  integration  limits  on  the  right  hand  side  of  Equation  3.74  is 

through  the  two  vertical  boundaries  on  the  left  and  right  hand  sides 

of  the  irregular  region  and  the  +  signs  represent  the  forces  that  are 

acting  to  or  away  from  the  finite  element  region.  All  in  all,  there 

are  a  maximum  of  four  different  forces  acting  at  the  vertical  boundaries 

thatare  determined  from  the  right  hand  side  of  Equation  3.74.  One  of 

these  forces  was  determined  by  Equation  3.69.  If  the  incident  energy  is 

from  the  left  then  there  are  three  forces  acting  at  the  vertical  boundaries. 

A  force  determined  by  Equation  3.69  due  to  the  incident  energy  at 

the  left  hand  boundary,  a  force  Pj^  due  to  reflection  of  energy  at  the 

left  boundary,  and  a  force  Pj  due  to  transmission  of  energy  at  the 

right  boundary.  Thus,  the  right  hand  side  of  Equation  3.74  is  integrated 

to  produce  the  above  three  forces,  for  an  arbitrary  fiV^.  Substitution 

of  Equations  3.6,  3.25a,  and  3.75  into  Equation  3.74,  the  first  term 

*T 

on  the  left  hand  side  is  integrated,  for  an  arbitrary  6V  ,  to  give 

/  N  D  B  dv  =  K  (3.76) 

V  — 

Matrix  is  an  8  by  8  stiffness  matrix  originating  from  four  nodal  points 
of  the  element,  each  node  having  a  horizontal  and  a  vertical  force.  To 
obtain  the  stiffness  matrix  for  all  elements  Equation  3.76  must  be  summed 
over  the  elements.  Thus, 

Kg  «  K  (3.77) 

Elements 

Substitution  of  the  contribution  from  Equations  3.45,  3.69,  and  3.77 
into  Equation  3.74  will  get 


I 

a: 
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(K  -  0,2  M)  V=  -  Pj_  .  (3.78) 

or 

(K  -  0,2m)  V  =  T  -  k ^  -  T  ^  -  L -  R(  \[  +  +  R^ 

(3.79) 

The  vector  V  represents  the  horizontal  and  vertical  displacements  of  the 
nodal  points  of  the  irregular  region.  Vfector  is  the  known  dis¬ 
placements  at  the  left  hand  boundary,  obtained  from  the  solution  of  the 
eigenvalue  problem.  The  unknown  vectors  the  reflected  displacements, 
and  vector  the  right  hand  boundary  displacement  can  be  incorporated 
into  vector  V  as  follows.  The  first  2N  elements  of  vector  V  are  equal 
to  vector  ^  and  if  the  vectors  are  substituted  for  the  2N 

A 

elements  of  vector  V,  then  the  unknown  vector  can  be  solved  for  by 

'il  ■  i  •  (S-SO) 

A  A 

The  last  2N  elements  of  vector  V  are  equal  to  vector  Vj..  Therefore, 
Equation  3.79  can  be  written  as 

(K  -  +  R  -  X  )  V  =  (R  -  L)  \[  .  (3.81 ) 

A 

The  unknown  vector  V  in  this  equation  can  be  solved  for  by  the  standard 
method  (28,  29). 


CHAPTER  IV 


ENERGY  OF  A  PROPAGATING  WAVE 
A.  Sumnary 

Elasticity  theory  assumes  the  medium  of  propagation  to  be  con¬ 
servative,  a  condition  in  which  no  energy  is  lost  or  changed  into  other 
forms.  In  reality  this  process  if  thermodynamically  reversible  only 
if  it  occurs  with  infinitisimal  speed,  so  that  the  thermodynamic 
equilibrium  is  established  in  the  body  at  every  instant.  Stressed  body 
however,  is  in  motion  and  has  finite  velocity.  The  body  is  not  in 
equilibrium  at  every  instant  and,  therefore,  processes  such  as  stress- 
induced  ordering  or  thermal  displacement  of  the  point  defect,  change 
in  dislocation  through  the  crystal  lattice,  and  viscous  sliding  at 
grain  boundaries  will  take  place  in  it  and  cause  it  to  return  to  equil¬ 
ibrium.  The  existence  of  these  processes  has  the  result  that  the  motion 
is  irreversible,  since  they  require  the  overcoming  of  a  barrier  to 
initiation.  At  the  end  of  each  cycle  of  stress  it  is  impossible  to 
reproduce  exactly  the  state  of  the  start;  hence,  there  will  be  hysteresis 
effects  and  energy  dissipation  or  attenuation  (Figure  1).  Because  of 
these  effects  the  amplitude  of  the  propagating  Rayleigh  wave  falls  off 
somewhat  more  rapidly  than  elastic  theory  predicts.  Although  attenuation 
is  extremely  important  in  practice,  it  must  still  be,  in  some  sense, 
small,  since  otherwise  waves  would  not  propagate  to  any  appreciable 
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distance.  For  large  attenuation  the  propagation  process  becomes  one 
of  the  diffusion  rather  than  propagation.  Since  waves 
propagate,  the  assumption  of  small  attenuation  keeps  the  theory  of 
elasticity  valid. 

The  dissipation  of  energy  arises  from  such  a  variety  of  mechanisms 
that  Is  Impossible  to  discuss  them  all  here.  The  most  effective  of 
all  mechanisms  fall  Into  the  class  of  viscoelasticity  In  which  the 
energy  of  the  wave  field  Is  taken  up  by  processes  such  as  defect  In 
mobility,  diffusion  rate,  phase  change,  etc.,  and  Is  restored 
partially  or  completely  to  the  field  with  a  delay.  As  a  result  of 
this  delay,  the  restored  energy  finds  itself  partially  out  of  phase 
with  the  original  wave  and  destructive  Interference  takes  place  (114). 

B.  Energy  Flux  Vector 

The  energy  flux  vector  E^(t)  at  time  t  1s  defined  as  the  rate 
at  which  energy  leaves  the  material  across  an  element  of  area  normal 
to  the  x^  -  axis,  measured  per  unit  area.  It  Is  given  by 

E.(t)  =  -  a.j(t)  Oj(t)  (4.1) 

where  Oj  is  the  velocity  vector  created  by  the  application  of  the  stress 
tensor  o^j,  the  minus  sign  Is  the  Indication  that  the  energy  Is 
leaving  the  material.  For  an  Isotropic,  elastic  medium,  the  direction 
of  the  energy  flow  Is  normal  to  the  wavefront.  Equation  4.1  fluctuates 
with  time  and  a  mean  energy  flux  vector  Is  defined  as  the  average  of 
E^(t)  over  a  cycle.  Thus, 
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in 

u  r  u) 


h-wT  E,(t)dt 


(4.2) 


The  energy  flux  vector,  E^,  contains  three  components.  Corres¬ 
ponding  to  the  directions  of  propagation,  i.e.,  in  vector  form. 


E^  .  E  .  c  .  U 


where  matrix  a  and  vector  U  are  given  by 


(4.3) 


®xx  °xy  °X2 

a  a  a 

yx  yy  yz 

COO 

2x  zy  2Z 


It  ^ '  It  ^  ‘  V  ■  V  ■ 


For  a  plane  wave  propagation  in  the  x  direction, 
vector  E  reduces  to  a  one  dimensional  energy  flux,  E,  and  is  given  by, 
setting  i  »  1,  and  j  =  1,3  in  Equation  4.1,  we  have 


(4.5) 


where  vectors  o  and  U  are  given  by 


®  *  ^%x*  ‘'xz^  •  It  ^  '  It  ^“x’  ^z^  EXP  (i(«t  -  Kx))  (4.6) 


the  stress  vector  represents  normal  and  shear  stresses  on  the  surface 
area  normal  to  the  direction  of  propagation.  The  energy  in  Equation  4.5 
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In  an  elastic  medium  Is  normal  to  the  wavefront,  I.e.,  the  energy  is  pro 
pagating  in  the  x  direction  normal  to  the  z  axis,  only  if  the  propa¬ 
gation  constant  K  is  real.  Let  us  analyze  the  state  of  vector  ic  in 
general.  From  Equation  4.4,  vector  K  in  a  plane  propagation  is 

ic  =  V  +  (4.7) 

A  A 

where  x  and  z  represent  direction  cosines  in  the  x  and  z  directions, 
respectively.  If  the  matrix  moduli  D  in  Chapter  III  contains  the 
effective  moduli  and  discussed  in  Chapter  II,  Section  C,  then 
vector  K  is  complex  and  is  given  by 


K  =  -  iKj  . 


(4.7a) 


Rewriting  the  vector  U  of  Equation  4.4  in  two  dimensions 

EXP  (i(u)t  -  ic  •  X)), 


U,(z) 


U 


where 


X  »  XX  +  zz, 

and  substituting  Equation  4.7  we  have 
U  ■  U(z)e  ^ 


(4.8) 


(4.9) 


(4.10) 


where  Kj  is  called  the  attenuation  vector.  To  include  attenuation  from 
the  vector  U(z)  it  is  better  to  write  Equation  4.10  as 


-k'.x 

U  -  |U(z)|  e  ^ 


gi(u)t  -  Kr*X) 


(4.11)  ' 
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The  dot  product  •  X  can  be  expanded  as 


•  X  =  eX  +  yZ  +  AR6(U(z))  (4.12) 


where  e  and  y  represent  the  components  of  Kj  in  x  and  z  directions, 
respectively,  and  ARG(U(z))  is  the  attenuation  factor  due  to  the  vector 
U(z).  The  existence  of  a  surface  wave  depends  on  whether  vector  Kj  is 
zero  or  not. 

The  attenuation  vector  Kj,  given  by  Equation  4.7a,  is  not,  in 
general,  perpendicular  to  the  propagation  vector  the  angle  e 
between  the  directions  of  the  two  vectors  is  between  zero  and  n/2. 


i.e.. 


0  <  e  <  n/2. 


(4.13) 


when  Q  =  0,  Kj  =  0,  and  therefore  there  exists  an  undamped  wave. 

When  0  =  n/2,  the  direction  of  attenuation  vector  Kj  is  perpendicular 
to  the  direction  of  propagation  vector  Kj^  and  the  wave  propagates  in 
the  direction  normal  to  the  surface  of  constant  phase  and  attenuates 
in  the  direction  normal  to  the  surface  of  constant  amplitude.  This 
occurs  when  6  =  y  =  0  in  Equation  4. '2  and  the  attenuation  is  due  to 
the  vector  U(z);  that  is,  the  amplitude  will  dampout  in  the  z 
direction.  This  is  the  classical  Rayleigh  wave. 
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The  surface  of  constant  phase  may  be  used  to  determine  the  phase 
velocity  of  the  wave  by  setting  the  progressive  part  of  Equation  4.11 
equal  to  a  constant,  say,  n/2. 


or 


gi(wt  -  K-X) 


»e"/2 

(4.14) 

-  in/2 

(4.15) 

Performing  the  dot  product  in  Equation  4.15  we  get, 

ut  -  KX  -  K, Z  =  -  in/2  . 

^  2 


(4.16) 


Taking  the  partial  derivative  of  x  and  z  with  respect  to  t  in  Equation 
4.16,  we  have, 

(4.17) 

Solving  for  and  ,  the  phase  velocities  in  the  horizontal  and 

O  w  V  V 

vertical  directions  are  given  by 


•  W 


•x  K, 


U)  y  ^  0) 

1/  » 


Z  K, 


(4.18) 


x  ■  z 
If  the  vector  K  is  complex  then 


Kx  =  Kx  ^  +  iy  , 


(4.19) 


and 


i*i(K  -is)  «K  o 


(4.20) 


or  approximately 


-  '»>• 


(4.21) 
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for  small  6  such  that  e 


2 


0. 


Similarly, 


=  V^d  -  iy).  (4.22) 

2 

for  small  y  such  that  y  0.  The  size  of  b  and  y  are  ordinarily  between 
.01  and  .001  of  the  and  of  the  K^,  respectively  so  that  the  above 
assumption  is  valid. 


C.  Energy  Flux  in  Terms  of  the  Wave' 
Number  and  the  Amplitude 


Let  us  rewrite  Equation  4.5  and  integrate  it  through  a  period  as 
in  Equation  4.2;  the  result  is 

E  =  /  “  E  dt  .  /  “  (|^  u)^.  ;  dt  (4.23) 

0  0 


Since  most  of  this  study  is  concerned  with  the  energy 
and  its  balance,  the  analysis  of  the  derivation  of  the  energy  relationship 
seems  essential.  Thus,  we  study  Equation  4.23  in  detail.  Consider  the 
dot  product  U  •  a  which  is  the  energy  density  normal  to  the  z  plane; 
the  energy  through  that  plane  is 

a  dz  =  Dy  e  dz  (4.24) 

0  0  — 


where 


X  0 
0  u 


(4.25) 


and 


i 
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Decomposing  the  matrix  in  Equation  4.25  as 


(4.25a) 


(4.26) 


and  substituting  this  and  Equation  4.25  into  Equation  4.24  gives 
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Substituting  Equations  3.38»  3.40  and  3.42  Into  Equation  4.27  we 


have 


/^U^-D„ed2  *  N  V  dz  -  IK  /V-N^  a  N  V  dz  (4.28) 

0  —  0  0 


where 


c-r  n 

L  u  0—*  L  0  v 


(4.29) 


From  Equation  3.42,  If  we  substitute 


1.  «  _L  iL 

3z  h/2  3v 


(4.30) 


and 


-  h  L  0  -1  0  1  J 


(4.31) 


Into  Equation  4.28  and  changing  the  variable  of  Integration  we  have 


h.i 


1.1 


/  U  -0  cdz  *  /  V'l/2 
0  —  -1 


l-n  0 
0  l-n 
1+n  0 

.  0  l+nj 


1  r-1  0  1  01 

C  — 

"  ^  L  0  -1  0  1  - 


0  1 


-  IK  /^V^l/2 


l-n  0 
0  l-n 
1+n  0 

0  1+n 


,  ri-r,  0  l*n  0  1  Jh 

L  0  l-n  0  l-HiJ  ^ 


dn 


(4.32) 


or 
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where 


and 


S  * 


h., 

/  U '•  DyCdz  » 

0  — 

T  R  V  dn  -  iK  /v’’’  i  h  s  V  dn 
-1  ^  “  -1  ®  “ 

(4.33) 

0 

-X(l-n) 

0  Xd-n)  ‘ 

R  » 

“w(l“n) 

0 

0 

-x(l+n) 

w(l-n)  0 

0  X{l+n) 

(4.34) 

-u(l+n) 

0 

u(l+n)  0 

‘(l-n)2(x+2v) 

0 

{1-t,^)(x+2p)  0 

- 

0 

(1-ti)^(X+2u) 

(l"n)^w 

0 

0  (1-n 

(l+n)^(x+2ii)  0 

')v 

(4.35) 

L 

0 

(l-n^)p 

0  (l+n)^uj 

Integration  of  Equation  4.33  yields 
h. 


/  u'-  D  edz  =  v'.  C.  V  +  iK  v''^-  A.  V  -  v'.  (C.  +  iKA.)  V  (4.36) 

A  ^  V  J  J 


for  the  layer.  Where 


"O  X  0  -X 

2(x+2p)  0  X+2p  0 

p  0  -p  0 

h. 

,  and  A.  *  ^ 
— 0  0 

0  2w  0  p 

0  X  0  -X 

X+2p  0  2(x+2p)  0 

V  0  -p  0 

_  0  p  0  2p_ 

(4.37) 


Sunming  up  Equation  4.36  over  all  layers  results  in 
h-  .  - 

^  /  0'*’.  Dycdz  *  XI  ^  ^ 


Layers 


j-1 
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Recall  that 


U  -  .  and  G  «  N  V 


(4.39) 


therefore , 


V  «  and  V  «  iuV 


(4.40) 


Substituting  Equations  4.38  and  4.40  into  Equation  4.23,  with  a 
realistic  definition,  the  average  energy  which  is  carried  across  the 
vertical  boundary,  by  the  mode  K,  is  given  by 


E  =  lif  /“  R(i<oV)'  •  (C  +  iKA)  R(V)dt 


when  matrices  A  and  £  are  real ,  and  by 
2il 

E  «  fr  /“  R(i«V)^  •  (R(C)  +  iKR(A))R(V)dt 


(4.41) 


(4.42) 


when  matrices  A  and  C  are  complex,  where  R  represents  the  real  part. 


Now  let 


V  «  (V^  +  iV2)e''“^  »  (V^  +  iV2)(cos  «t  +  i  sin  ut)  (4.43) 


Then  substituting 


R(V)  ■  cos  ut  -  V2  sin  ii)t, 

A  A  A 

R(iuV)  ■  -  ta(V^  sin  ut  +  V2  cos  lot) 


(4.44) 


R(iwV)^*  R(V)  =  -  o)Vj  •  sin  wt  cos  ut  -  vj  •  V2  cos^  ut  + 


wVj  .  V2  sin^  ut  +  vj  •  Vg  sin  ut  cos  ut  (4.45) 
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2i  ^  2n 

Ijj-  /“  sin  »t  cos  wt  dt  «  0,  |j-  /“  sin^  wtdt  *  fjf  /“  cos^  wtdt  «  1/2 

(4.46) 

Into  Equation  4.41  yields, 

E  »  I  (vj  •  (C  +  m)V2  '  •  (£:+  i»CA)V^),  (4.47) 

where  and  V2  are  the  real  and  imaginary  of  the  vector  V. 

Equation  4.47  expresses  the  energy  carried  by  the  mode  K.  If  the 
incident  energy  of  the  Rayleigh  wave  across  the  vertical  boundary  is  in 
the  fundamental  mode  K,  given  by  Equation  4.47,  then  the  reflected  and 
transmitted  energies  must  be  calculated  through  the  modes  which  are 
produced  by  the  incident  energy.  Therefore,  if  L  is  the  number  of 
different  modes  produced  by  the  incident  energy  then  the  reflected  and 
transmitted  energies  are  given  by 

S=  I 

and 

Et  •  j,  ?  <*u-  <£  *  '£  *  >  <<•«> 

respectively. 

0.  Attenuation 

If  we  substitute  the  effective  moduli  of  Equation  2.17  into  the 
energy  relation  Equation  4.47  it  produces  an  additional  term  called 
the  power  loss  or  dissipated  power  in  the  system.  This  power  loss  is 
due  to  the  viscoelastic  part  of  Equation  2.17,  and  is  given  by 
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Pp  *  f  (vj-  (IM(C)  +  1KIM(A))V2  -  V2-  (IM(p  +  iKIM(A))V2),  (4.48) 


where  IM  represents  the  Imaginary  part  or  the  viscoelastic  part  of  the 
matrices  C  and  A.  The  ratio  of  the  average  power  loss.  Equation  4.48, 
to  the  average  rate  of  energy  transport.  Equation  4.47,  is  called  the 
attenuation  and  is  denoted  by  a  dimensionless  parameter  2nQ’^;  thus. 


p 

2jjq-1  _  average  dissipated  power  ^  _D 
^  ~  rate  of  energy  transport  E  ‘ 

If  the  complex  wavenumber 


(4.49) 


K  =  Kp  +  iKj  (4.50) 

is  given,  or  for  a  given  period,  is  found  from  the  solution  of  the  eigen¬ 
value  problem  for  the  layered  region,  then  the  attenuation  is  given  by 


Equation  4.49  is  general;  it  can  be  used  for  spatial  attenuation,  i.e., 
attenuation  arising  from  discontinuities  along  the  path  of  the  wave, 
or  for  temporal  attenuation,  i.e.,  attenuation  arising  from  time, 
temperature  or  other  viscoelastic  behavior  of  the  medium. 


E.  Change  of  Wavenumber  at  the  Vertical  Boundary 

The  existence  of  reflected  energy  at  the  vertical  boundary 
suggests  that  the  phase  velocity  at  the  left  and  the  right  hand  sides 
of  the  boundary  are  different.  Since  the  wave  is  harmonic  then  the 
change  in  velocity  is  due  to  the  change  in  wavenumber.  This  change 
must  be  calculated  by  means  other  than  from  Equation  3.81  whose  solution 
gives  the  displacements  at  the  vertical  boundary. 
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Au1d  (7,  8)  has  used  a  perturbation  theory  and  has  derived  a 
relation  to  calculate  the  change  in  the  wavenumber  due  to  various 
mechanisms.  Auld's  relation  may  be  revised  to  be  used  in  seismic 
waves;  the  form  suitable  here  may  be  expressed  as 


4K  =  u 


/  AD  c  -  U)  dz 

z 


4E 


where 


X+2v 


D  * 


X 


0 


X  0‘ 
X+2p  0 
0 


(4.52) 


(4.53) 


3X  3Z  / 


(4.54) 


G  -  G(z)e^^“^  ■ 


(4.55) 


E  is  given  by  Equation  4.49  and  t  denotes  the  contrast  between  the  two 
different  media.  Use  of  Equations  3.33  and  3.34  reduces  the  numerator 
of  Equation  4.52  to 


1“* 

O 

_ 1 

1 — 
a 

3X 

3X 

»  If 

G)  •  A  0 

»  If 

3  3 

i_  1- 

1 — 

CP| 

X  1 

_3Z  3X_ 

0  -  Apw^  G)  dz  .  (4.56) 


Decomposition  of  the  operating  martix  Into  two  matrices,  depending  on 
X  or  z  only,  gives 


50 


— » 

r 

o 

_J 

1 

r 

o 

~  “■ 

0  0 

3X 

ax 

0  ^ 

S 

0  0 

+ 

0  i- 

3Z 

3Z 

3_  ^ 

0  — 

1.  0 

_  3Z  3X  _ 

L  3X  J 

Laz  J 

(4.57) 


and  operating  on  the  matrix  with  the  partial  derivative  with  respect  to 
X  on  the  corresponding  vectors  results  In 


L  *  /((vz  -  iKI)U)^*AD  (TZ  -  iKI)U  -  pu)^U^‘  U)dz,  (4.58) 

z 

where  matrix  vz  Is  given  by  the  second  term  of  Equation  4.57,  and 


I 


1  0 
0  0 
Lo  Ij 


Again,  by  using 

U  =  N  V 
.'l-n 


0  1+n 

l-n  0 


0 

1+n 


,  -  1  <  n  <  1 


(4.59) 


(4.60) 


Z  3  1  3  -  . 

'’•-h72“*  3l'  Oiz  <h 


fz-^r  "  1  ’ 

"  L  0  -1  0  iJ 


(4.61) 


(4.62) 


Equation  4.58  becomes  (for  the  moment  we  set  AD  »  D,  Ap  *  p) 
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L  « 


0 

0 

■1 


— 1 

o 

o 

o 

o 

_ 1 

1-n  0 

1+n 

0 

0-101 

- 

IK/2 

0  0 

0 

0 

-10  10 

0  1-n 

0 

1+n 

— 

0  0  0 

1-n 

0  1+n 

0 

1  0  1 

-  IK/2 

0 

0  0 

0 

)v 

0  1  0 

0 

1-n  0 

1+n 

■(1-n)^ 

0 

l-n^ 

0 

0 

(1-n)^ 

0 

l-n2 

\  - 

5 

0 

Vh/2dn. 

1-n^ 

0 

(1+n) 

^  0 

/ 

0 

1-n^ 

0 

(1+Tl)^ 

2  .T 

U  II I 

p  -g-  V  • 


After  performing  the  transpose  and  corresponding  multiplications, 
Equation  4.63  becomes 


L-  /  (^”^.(12 
-1 


u  0  “u  0 

0  a+2w  0  -(A+2u) 

-u  0  w  0 

0  -(x+2u)  0  X+Zv 


-  i 


2F 


0 

-A(l-n) 

0 

X(l-Tl) 


■u(l-n) 

0 

u(l-n) 

0 


•M(l+n) 

0 

u(l+n) 


(4.63) 


-A(l+n) 

0 

x(l+n) 


0 
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0 

-x(1-n) 

0 

x(l-n) 

K 

-uC^-n) 

0 

p(l-n) 

0 

2h 

0 

-w(l+r\) 

-x(l+n) 

0 

0 

u(1+ri) 

x(l+ti) 

0 

(x+2u)(1-n)‘ 

0 


0  {x+2y){1-n^)  0 

I  \  2  jh  / « 


.2  0  u(l-n)^  0  w(l-n^)  v  . 

(x+2u)(1-t)  )  0  (x+2y)(l+n)^  0  ^ 


p(1"n^) 


p(1+n)' 


2  -T  0 

-  P  F  V^.  , 


(1-n)  0 


0  (1-n)^  0 


0  (1+n)^  0 

l-n^  0  (l+n)‘ 


Vh/2dn  (4.64) 


Now,  the  integrations  of 
1  1 


/  (l-n)eln  *  /  (1+Ti)dri  *  2,  /  (l-n)^dr)  *  /  (l+n)^dTi  =  4  ,  (4.65) 

-1  -1  -1  -1 


/  (l“n  )<ln  “  7 
-1  ^ 


reduces  Equation  4.64  into 


L  »  V^*  (V2  G.  +  iKF  +  K^A.  -  w^M.)  V 

W  V  V 


(4.66) 


where  matrices  Gj,  Aj,  and  Mj  are  given  by  Equations  3.48  -  3.50  and 


-(x-u) 

0 

(x+p) 

0 

X+p 

0 

“(x+p) 

0 

X-p 

-(X-w) 


(4.67) 


Since  j  represents  a  typical  layer,  the  sum  of  Equation  4.66  over  all 
layers,  when  substituted  into  Equation  4.52  gives,  considering  now  aD 
and  aP, 


6K  =  (1) 


V^-(l/2  AG  +  iKA  F  +  K^A  A  -  Jt  M)V 


(4.68) 


When  calculated,  this  AK  may  be  added  to  the  wave  number  of  the  left 
hand  side  of  the  boundary;  the  phase  velocity  of  the  right  hand  side  of 
the  boundary  can  be  calculated  by  the  expression. 


V  s  *** 

\  K+6K  • 


(4.69) 


For  different  periods  the  dispersion  curve  of  the  desired  medium  can  be 
computed. 

Equation  4.68  may  also  be  used  for  the  determination  of  attenuation. 
For  this  purpose  if 


«K  =  +  iAKj 

R 


(4.70) 


AKj  *  u 


(1/2  aGj  +  iK  a  Fj  +  a  Aj)V 


and  the  attenuation  is  given  by 


(4.71) 


where  I  stands  for  the  Imaginary  part.  Since  attenuations  due  to  various 
mechanisms  are  additive.  Equation  4.72  may  be  added  to  the  attenuation 
obtained  due  to  viscoelastic  behavior. 
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CHAPTER  V 


APPLICATION  OF  THE  METHOD 
A.  A  Model 

To  apply  a  model  to  the  theories  discussed  in  the  last  chapters, 
we  chose  Herrin's  models  (49)  of  the  Canadian  Shield  and  Basin  and  Range 
provinces  for  the  elastic  part  of  tht  medium.  The  Basin  and  Range 
province  is  used  in  the  left  layered  region  and  the  Canadian  Shield  is 
used  in  the  right  layered  region;  the  order  is  reversible.  For  the 
irregular  region  we  chose  different  combinations  of  the  Basin  and  Range 
province  and  the  Canadian  Shield.  The  data  given  in  Herrin's  model  are 
thicknesses,  densities,  compressional  velocities,  and  shear  velocities 
of  the  two  layered  media.  To  account  for  the  anelasticity  in  the  model 
a  complex  part  for  the  compressional  velocities  and  shear  velocities 
may  be  assumed.  Since  the  elastic  moduli  are  obtained  by  using  the 
expressions 


X  •  «(Vp  -  2v/) 

(5.1) 

u  ■  , 

(5.2) 

where,  \  and  u  are  the  elastic  moduli  of  the  medium,  and  Vp  and  are 
the  real  parts  of  compressional  and  shear  velocities,  respectively, 
then  the  effective  moduli,  or  complex  moduli,  as  given  in  Equation  2.17, 
may  be  obtained  by 
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p(Vp^  - 

(5.3) 

(5.4) 

for  complex  velocities.  For  this  reason,  in  our  model,  any  viscoelastic 
behavior  of  the  medium  may  be  incorporated  in  the  imaginary  parts  of  the 
compressional  and  shear  velocities.  Thus,  for  the  viscoelastic  part 
of  the  medium  in  this  paper,  three  different  sets  of  data  are 
used:  A  small  enough  (.1  percent  of  the  real  part)  imaginary  parts  for 
compressional  and  shear  velocities;  a  starting  constant  shear  quality 
factor,  Qg,  for  the  Basin  and  Range  province  and  a  starting  constant  Qg 
for  the  Canadian  Shield;  and  finally  the  starting  Qg  model  of  Lee  and 
Solomon  (69a),  whose  Qg  of  Western  U.  S.  and  East-Central  U.  S.  are 
used  in  Basin  and  Range  provinces  and  Canadian  Shield,  respectively. 

When  Qg  is  given,  the  imaginary  parts  of  Vp  and  are  determined 


by  the  following  expressions. 

From 

REAL  (Vj) 

^6  '  2  TMAG  (Vj) 

(5.5) 

Vp  =  REAL  (Vp)  +  i 

IMAG  (Vp), 

(5.6) 

Vg  =  REAL  (Vj)  +  i 

IMAG  (Vj), 

(5.7) 

REAL  (Vj 

(5.8) 

IMAG  (Vj)  = 

And  from  {69a), 
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we  have 


IMAG(Vp) 


REAL  (Vp) 


(5.9) 


(5.10) 


Both  regions  of  the  model  are  assumed  each  to  contain  15  layers 
with  a  total  depth  of  250  kilometers;  this  depth  includes  the  crust  and 
upper  mantle  so  that  any  partial  melting  zone  in  the  upper  mantle  may  be 
studied.  The  irregular  region  I  is  divided  into  15  by  15  network  of 
elements  whose  total  number  of  elements  and  total  number  of  nodal  points 
of  the  elements  are  N  and  N(N  +  1),  respectively,  where  N  is  the  number 
of  layers.  The  slightly  adjusted  Herrin  model  for  the  above  geometry 
is  given  in  Tables  1  and  2.  The  integer  N  in  Table  2  represents  the 
number  of  elements,  in  the  finite  element  region,  that  have  the  same 
property.  The  numbers  assigned  to  the  integer  N  in  Tables  2  and  2A  are 
completely  arbitrary. 

B.  Real  Data 

A  set  of  periods  from  9  to  60  seconds  along  with  the  data  of 
Table  1,  may  be  used  as  an  input  for  the  solution  of  the  eigenvalue 
problem.  Equation  3.52,  in  the  layered  regions.  The  fundamental  solution 
of  the  eigenvalue  problem  corresponding  to  each  period  is  the  largest 
eigenvalue  or  the  largest  wavenumber  selected  from  the  2N  eigenvalues. 

The  phase  velocity  corresponding  to  each  period  is  given  by 


PHASE  VELOCITY  = 


OMEGA 

WAVENUMBER  ’ 


(5.11) 


for  the  given  period  and  corresponding  frequency. 


The  computer  program  RALEE  used  to  calculate  the  unknown  of  the 
medium  is  the  revised  form  of  the  program  RALEEM  written  by  Drake  (29) 


in  1972.  To  test  the  program  RALEE  we  used  Lysmer  and  Drake's  (80) 
data  of  the  Sierra  Nevada  and  calculated  the  phase  velocities  of  Rayleigh 
waves  for  a  range  of  periods  of  6  to  10  seconds  and  the  horizontal  and 
vertical  amplitudes  for  a  period  of  6  seconds.  The  agreement  is  excel¬ 
lent  for  the  phase  velocities  and  good  for  the  amplitudes;  but  the 
small  difference  in  the  amplitudes  may  be  due  to  scale  arbitrariness. 

The  phase  velocities  are  given  in  Table  16  and  the  amplitudes  are  shown 
in  Figure  16. 

The  wavenumbers  and  the  phase  velocities  calculated  for  the  two 
layered  systems  of  the  Basin  and  Range  province  and  the  Canadian  shields 
with  the  data  of  Table  1  and  the  periods  of  9  to  60  seconds  are  also 
given  in  Tables  3  and  4,  and  the  corresponding  dispersion  curves  are 
shown  in  Figures  7  and  8. 

To  test  the  effect  of  layering  we  reduced  the  number  of  layers 
from  15  to  3  and  assumed  two  layers  for  the  crust  and  a  deep  layer  for 
the  upper  mantle.  The  computed  Rayleigh  wave  phase  velocities  are  some¬ 
what  higher  than  those  of  the  15  layered  case.  The  difference  may  be 
due  somewhat  to  the  averaging  data  of  15  layers  into  3  layers  but  it 
is  mostly  due  to  the  large  dimensions  of  the  elements  as  compared  to 
the  wavelength  of  the  shear  wave  (104). 

To  test  the  computer  program  RALEE  for  the  irregular  zone  I  in 
the  real  case,  we  first  assumed  the  three  zones  L,  I  and  R  of  Figure  5 
to  have  the  same  properties  as  the  Basin  and  Range  provinces  and  calcu¬ 
lated  the  reflected  and  transmitted  energies.  For  the  period  range  of 
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10  to  80  seconds  the  sum  of  the  reflected  and  transmitted  energies  is 
very  close  to  incident  energy.  These  data  are  shown  in  Table  5. 

Second,  we  set  the  boundary  between  the  Basin  and  Range  province  and 
the  Canadian  Shield  at  the  right  hand  side  of  the  tone  L,  i.e.,  we  let 
the  zone  L  be  the  Basin  and  Range  province  and  zones  I  and  R  be  the 
Canadian  Shield.  Third,  we  set  the  boundary  at  the  left  hand  side  of  the 
zone  R,  i.e.,  we  let  the  zones  L  and  I  be  the  Basin  and  Range  province 
and  zone  R  be  the  Canadian  Shield.  Fourth,  we  set  the  boundary  at  the 
middle  of  zone  I,  i.e.,  we  let  the  zone  L  and  one-half  of  the  zone  I  be 
the  Basin  and  Range  province  and  one-half  of  the  zone  I  and  zone  R  be 
the  Canadian  shield.  Finally,  we  let  the  zones  L  and  R  be  the  Basin 
and  Range  province  and  the  Canadian  shield,  respectively,  but  gradually 
changed  the  parameters  of  the  zone  I  from  the  parameters  of  the  Basin 
and  Range  province  into  the  parameters  of  the  Canadian  shield.  The  sum 
of  the  reflected  and  transmitted  energies  determined  for  each  of  the 
above  models  in  the  real  case,  was  found  to  be  almost  equal  to  the 
incident  energy.  The  transmitted  energy  in  all  of  the  above  models, 
in  the  real  case,  is  more  than  99%  except  for  a  few  period  of  18  to  25, 
in  the  second  case,  in  which  for  each  period  the  reflected  energy  is 
the  highest,  i.e.,  from  .05%  to  2.2%.  On  the  other  hand,  in  the  third 
case,  the  reflected  energy  is  the  smallest,  i.e.,  from  .007%  to  .06%. 
Therefore,  the  decrease  in  reflected  energy  is  from  left  to  right  with 
respect  to  the  position  of  the  vertical  boundary;  nevertheless,  the 
process  is  reversible.  The  values  for  the  above  cases  are  given  in 
Tables  6,  6A,  6B,  and  6C,  respectively. 


For  the  determination  of  the  phase  velocity  in  the  irregular 
region,  in  the  real  case,  the  boundary  displacements,  determined  by 


Equation  3.81,  are  used  in  Equation  4.68;  the  left  boundary  displace¬ 
ments  are  replaced  by  the  vector  and  the  right  boundary  displace¬ 
ments  are  replaced  by  the  vector  V  in  Equation  4.68.  Substitution  of 
the  resulting  into  Equation  4.69  determines  the  fundamental  mode 
phase  velocity.  Repetition  of  the  use  of  Equations  4.68  and  4.69, 
corresponding  to  the  number  of  distinct  periods  gives  the  phase 
velocities  of  the  irregular  region  for  each  case  mentioned  in  the 
previous  paragraph.  The  dispersion  curves  for  each  case  are  shown  in 
Figures  9,  9A,  98,  9C,  and  9D,  respectively. 

C.  Complex  Data 

The  geometry  of  the  model  remains  the  same  in  the  viscoelastic 

case.  Therefore,  the  real  parts  of  the  complex  V  and  V  velocities  do 

P  5 

not  change  in  the  complex  case.  Since  thicknesses,  densities,  and 
periods  are  assumed  to  be  real  in  the  elastic  and  viscoelastic  medium, 
they  will  not  be  repeated  again  in  this  section.  Thus,  for  the  complex 
case  of  the  model  we  assigned  three  different  sets  of  data  to  the 
imaginary  parts  of  the  complex  Vp  and  velocities.  First,  the  imag¬ 
inary  parts  of  the  V  and  V  were  assumed  to  be  .1  percent  of  their 
P  ^ 

corresponding  real  parts.  Second,  a  constant  shear  quality  factor  for 
each  zone  of  L  and  R  was  assumed  and  the  corresponding  compressional 
quality  factor  and  the  imaginary  parts  of  Vp  and  were  calculated. 
Third,  we  used  the  Q.  model  of  Lee  and  Solomon  (69a)  and  calculated 

P 

the  corresponding  compressional  quality  factor  and  the  imaginary 
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parts  of  and  V^.  The  imaginary  parts  of  and  V,  for  the  three 
pi  p  s 

cases  are  shown  in  Tables  7,  8,  and  9. 

For  the  same  number  of  periods  that  was  used  in  the  real  case, 
the  values,  given  in  Tables  7,  8,  and  9  are  used  as  the  input  for  the 
solution  of  complex  eigenvalue  problem.  Equation  3.52,  in  the  layered 
regions.  The  fundamental  solution  of  the  complex  eigenvalue  problem 
corresponding  to  each  period  is  the  largest  eigenvalue  with  a  negative 
imaginary  part,  selected  from  2N  eigenvalues.  The  phase  velocity  and 
the  attenuation  corresponding  to  each  period  are  given  by 


PHASE  VELOCITY  = 


OMEGA _ 

REAL  (WAVENUMBER) 


.  2  X  IMAG  (WAVENUMBER) 
AHENUATION  *  imiumm  ’ 


(5.11a) 


(5.12a) 


respectively.  The  real  and  imaginary  parts  of  wavenumbers,  the  phase 
velocities  and  the  attenuations,  corresponding  to  Tables  7,  8,  and  9 
are  given  in  Tables  10  -  15  and  the  corresponding  attenuation  curves 
and  the  dispersion  curves  are  shown  in  Figures  10  -  15. 

For  the  complex  part  of  the  irregular  region  I  we  used  the 
Imaginary  parts  for  V„  and  V,  velocities,  obtained  from  the  Q*  of  Lee 

p  S  p 

and  Solomon  (69a),  in  all  cases  described  for  the  real  case  in  section  B. 
In  each  case  the  sum  of  the  reflected  and  transmitted  energies  is 
considerably  less  than  the  incident  energy,  but  the  difference  is 
remarkably  in  accordance  with  the  given  Q.  value  with  respect  to  the 

P 

depth  where  the  energy  is  computed.  The  attenuated  percentage  of 
reflected  and  transmitted  energies  are  given  in  Tables  5,  6,  6A,  6B, 
and  6C. 


For  consistency  of  the  result,  we  asked  H.  Mack  (80a)  to  input 


our  complex  data  In  his  computer  program  VERES.  He  calculated  the 
transmitted  Rayleigh  wave  energy  for  the  case  where  the  medium  is 
considered  to  have  the  parameters  of  Basin  and  Range  province.  His 
result  agrees  well  with  higher  periods,  fluctuates  somewhat  with  the 
intermediate  periods,  and  agrees  with  periods  of  13  -  15  seconds  as 
shown  in  Table  5.  For  periods  lower  than  13  seconds  his  program  is 
not  applicable. 

The  determination  of  the  phase  velocity  of  the  irregular  region, 
in  the  complex  case,  follows  the  same  procedure  of  the  last  subsection; 
however,  since  the  change  in  the  phase  velocities  between  the  real  and 
the  complex  cases  are  very  small  the  dispersion  curves  would  look  the 
same  as  the  ones  i>i  the  real  case.  Thus,  the  repetition  is  not  necessary. 

Finally,  for  a  check  on  the  result  we  used  the  data  of  Lee  and 
Solomon  (69a,  pages  77  and  83)  in  our  15  layered  systems  and  computed 
the  phase  velocities  and  the  attenuations  for  East-Central  and 
Western  North  America.  The  result  is  compared  well  with  their  result 
which  is  obtained  by  inversion.  The  phase  velocities  and  the  attenuations 
are  given  in  Tables  17  and  18  and  the  corresponding  dispersion  and 
attenuation  curves  are  shown  in  Figures  17,  18,  19,  and  20. 


CHAPTER  VI 


SUMMARY  AND  CONCLUSION 

The  anomalous  decay  of  surface  wave  amplitude  may  be  due  to 
viscoelasticity  and  geometrical  discontinuities  of  the  medium.  In 
this  paper  the  path  of  a  Rayleigh  wave  was  assumed  to  be  through  the 
Basin  and  Range  provinces  and  the  Canadian  Shield.  The  two  regions 
were  treated  as  elastic  and  inelastic;  and  the  surface  of  intersection 
between  them  was  considered  to  be  an  approximation  to  a  plane  vertical 
boundary. 

With  the  plane  vertical  boundary  removed*  the  two  regions,  as  a 
layered  system,  were  studies  as  separate  elastic  and  viscoelastic  media. 
The  Rayleigh  wave  phase  velocities,  amplitudes,  and  the  attenuations 
of  the  two  regions,  with  the  technique  of  finite  element,  were  computed 
for  a  range  of  periods.  The  change  in  phase  velocities,  in  each  region 
due  to  viscoelasticity  in  three  different  cases  of  data,  were  from  .01 
percent  to  0.5  percent  of  elastic  phase  velocities,  an  almost  negligible 
amount.  Consequently,  the  dispersion  curves,  with  the  range  of  periods 
of  9  -  60  seconds,  are  almost  similar  in  elastic  and  viscoelastic  cases 
in  all  three  zones.  On  the  other  hand,  the  amplitudes  and  attenuations, 
in  three  cases  of  different  viscoelastic  parameters  in  the  layered 
system,  are  in  accordance  with  the  input  data  and  agree  well  with  the 
very  limited  published  data.  Though  small,  the  difference  in  amplitude 
and  attenuation  due  to  viscoelasticity  appear  to  be  significant  and  is 
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1n  accordance  with  the  amount  of  imaginary  parts  of  and  velocities 
with  respect  to  the  depth  of  the  layers.  The  result  of  the  use  of  the 
model  of  Lee  and  Solomon  (69a)  also  confirms  the  above  result,  in 
addition  to  being  in  excellent  agreement  with  their  result. 

The  reflected  and  transmitted  energies  through  the  vertical 
boundary  are  functions  of  position  of  the  vertical  boundary  with  respect 
to  the  two  media  of  L  and  R,  in  the  real  case.  They  are  in  addition 
to  the  above  condition,  functions  of  attenuation  or  the  Q.  as  well, 

P 

in  the  viscoelastic  case.  The  reflection  of  energy,  although  small  in 
all  cases,  decreases  as  the  vertical  boundary,  in  our  model,  moves  from 
left  to  right;  but  the  process  is  completely  reversible.  The  highest 
reflection  of  energy  (.04%  to  2.2%)  occurs  where  the  vertical  boundary 
is  at  the  right  hand  side  of  the  zone  L  so  that  the  zones  I  and  R  both 
have  the  parameters  of  the  Canadian  shield  and  the  zone  L  has  the  Basin 
and  Range  parameters.  The  one  position  of  the  vertical  boundary  where 
the  reflected  energy  has  the  intermediate  value  (.04%  to  .5%)  is  in 
the  middle  of  the  zone  I;  the  reflected  energy  in  the  real  and  complex 
cases  are  almost  equal;  the  sum  of  the  reflected  and  transmitted 
energies,  in  the  real  case  are  almost  exactly  equal  to  the  incident 
energy  and  the  transmitted  energy,  in  the  complex  case  is  the  highest. 

Due  to  the  fact  that  the  reflected  energy,  in  all  of  the  models 
tested,  is  very  small,  it  may  be  concluded  that  for  a  normal  incidence 
of  Rayleigh  wave  at  a  vertical  boundary  between  the  Continental  provinces 
the  mode  conversion  is  too  small  to  measure  in  the  real  crust.  Due 

to  the  computer  core  limitation  the  number  of  sublayers  and 
the  size  of  the  elements  may  be  somewhat  large  compared  to 
the  wavelength  of  the  shear  wave. 
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TABLE  1 

HERRIN'S  MODEL  OF  BASIN  AND  RANGE  PROVINCES 
AND  CANADIAN  SHIELD,  SLIGHTLY  ADJUSTED  TO 
FIT  THE  FIFTEEN  LAYERED  MODEL 


Basin  and  Ranae  j 

Canadian  Shield 

t _ 

Depth 

(KM) 

Density 

(G/Cm3) 

Vp 

(KM/SEC) 

(KM^SEC) 

Depth 

(KM) 

Density 

(G/CM3) 

Vp  Vs 

(KM/SEC)  (KM/SEC) 

0-14 

2.640 

6.000 

3.540 

0-14 

2.610 

5.800 

3.470 

14-19 

2.640 

6.000 

3.560 

14-19 

2.620 

6.700 

3.740 

19-20 

2.640 

6.700 

3.600 

19-20 

2.620 

6.700 

3.740 

20-28 

2.640 

6.695 

3.595 

20-28 

2.620 

6.700 

3.740 

28-30 

2.640 

6.650 

3.550 

28-30 

2.620 

6.700 

3.740 

30-34 

2.900 

7.900 

4.430 

30-34 

2.620 

7.000 

3.820 

34-50 

2.900 

7.900 

4.425 

34-50 

2.900 

8.050 

4.620 

50-60 

3.420 

7.899 

4.420 

50-60 

2.900 

8.050 

4.620 

60-80 

3.400 

7.880 

4.390 

60-80 

2.900 

8.050 

4.620 

80-90 

3.400 

7.860 

4.355 

80-90 

2.920 

8.100 

4.540 

90-100 

3.392 

7.820 

4.320 

90-100 

3.400 

8.400 

4.660 

100-120 

3.382 

7.678 

4.182 

100-120 

3.500 

8.500 

4.500 

120-140 

3.380 

7.730 

4.220 

120-140 

3.500 

8.500 

4.500 

140-150 

3.400 

7.800 

4.270 

140-150 

3.500 

8.500 

4.500 

150-200 

3.400 

7.890 

4.430  ! 

150-200 

3.500 

8.500 

4.500 
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TABLE  2a 

DATA  FROM  TABLE  2,  WITH  LEE  AND  SOLOMON'S  DATA,  FOR 
THE  FINITE  ELEMENT  REGION  IN  COMPLEX  CASE 


Vp  (KM/SEC) 

Vj  (KM/SEC) 

REAL 

(Vp) 

IMAG  (Yp) 
(xlO-^*) 

N 

IMAG  (Vc) 
(xlO-2) 

N 

6.00 

0.40 

15 

3.54 

6.00 

15 

5.80 

0.14 

1 

3.47 

0.17 

1 

6.00 

0.40 

14 

3.56 

6.00 

14 

6.70 

0.14 

2 

3.74 

0.18 

2 

6.70 

0.40 

13 

3.60 

6.00 

13 

6.70 

0.14 

4 

3.74 

0.18 

2 

6.69 

0.20 

3.59 

0.23 

13 

6.70 

0.14 

6  ' 

3.74 

0.18 

4 

6.66 

0.20 

9  i 

3.55 

0.23 

11 

6.70 

0.14 

8 

3.74 

0.18 

6 

7.90 

0.20 

7 

4.43 

0.28 

9 

7.00 

0.14 

10 

3.82 

0.18 

8 

7.90 

0.20 

5 

4.42 

0.28 

7 

8.05 

0.17 

12 

4.62 

0.23 

10 

7.90 

0.20 

3 

4.42 

0.28 

15 

8.05 

0.17 

14 

4.62 

0.23 

12 

7.88 

0.20 

1 

4.39 

0.28 

3 

8.05 

0.17 

11 

4.62 

0.28 

14 

7.86 

9. 00 

4 

4.36 

13.00 

1 

8.10 

0.17 

13 

4.62 

0.23 

11 

7.82 

9.00 

2 

4.32 

13.00 

4 

8.40 

0.16 

13 

4.62 

0.23 

13 

7.68 

9.00 

2 

4.18 

12.00 

2 

8.50 

0.16 

13 

4.62 

0.23 

13 

7.73 

9.00 

2 

4.22 

12.00 

2 

8.50 

3.60 

13 

4.54 

0.23 

13 

7.80 

9.00 

2 

4.27 

13.00 

2 

8.13 

3.60 

13 

4.66 

5.40 

13 

7.85 

5.00 

2 

4.30 

1 _ 

7.00 

2 
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TABLE  3 

FUNDAMENTAL,  AND  FIRST  AND  SECOND  HIGHER  MODE  EIGENVALUES 
AND  CORRESPONDING  PHASE  VELOCITIES  OF  THE 
LEFT  layers  (REAL  CASE) 


Wavenumbers  (KM" 

M  * 

Phase  Velocity  (KM/SEC) 

Funda¬ 

mental 

Mode 

1st 

Higher 

Mode 

2nd 

Higher 

Mode 

Funda¬ 

mental 

Mode 

1st 

Higher 

Mode 

2nd 

Higher 

Mode 

0.2009 

0.1618 

0.1589 

3.475 

4.315 

4.394 

0.1806 

0.1453 

0.1420 

3.478 

4.323 

4.423 

0.1497 

0.1207 

0.1168 

3.496 

4.337 

4.481 

0.1179 

0. 0960 

— 

3.551 

4,363 

— 

0, 0962 

0.0792 

— 

3.627 

4.402 

— 

0.0853 

~ 

— 

3.681 

— 

— 

0.0765 

~ 

— 

3.732 

~ 

— 

0. 0662 

~ 

— 

3.796 

~ 

— 

0.0541 

— 

— 

3.865 

— 

— 

0.0460 

— 

— 

3.902 

— 

— 

0.0400 

— 

— 

3.924 

— 

“ 

0.0354 

— 

— 

3.940 

~ 

— 

0.0318 

— 

— 

3.955 

— 

— 

0.0262 

— 

— 

3.997 

— 

— 

0.0221 

“ 

— 

4.066 

— 

— 
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TABLE  4 

RINDA^^EMTAL,  AND  FIRST  AND  SECOND  HIGHER  MODE  EIGENVALUES 
AND  CORRESPONDING  PHASE  VELOCITIES  OF  THE 


RIGHT  LAYERS 

(REAL  CASE) 

Wavenumbers  (KM‘ 

•') 

Phase  Velocity  (KM/SEC) 

Funda¬ 

1st 

2nd 

Funda¬ 

1st 

2nd 

mental 

Higher 

Higher 

mental 

Higher 

Higher 

Mode 

Mode 

Mode 

Mode 

Mode 

Mode 

0.2002 

0.1545 

0.1533 

3.487 

4.522 

4.553 

0.1799 

0.1368 

0.1321 

3.493 

4.537 

4.590 

0.1490 

0.1150 

— 

3.514 

4.552 

— 

0.1174 

0.0915 

“ 

3.569 

4.577 

— 

0.0958 

0.0757 

— 

3.645 

4.614 

— 

0.0849 

— 

— 

3.700 

— 

— 

0.0761 

— 

3.755 

— 

— 

0.0656 

— 

— 

3.830 

— 

— 

0. 0533 

— 

3.928 

— 

— 

0. 044  9 

“ 

~ 

3.996 

— 

— 

0.0388 

— 

— 

4.045 

— 

— 

0.0342 

~ 

— 

4.084 

— 

— 

0.0305 

4.120 

— 

— 

0.0250 

— 

— 

4.198 

— 

— 

0.0208 

— 

— 

4.308 

— 

— 

70 


IT 


TABLE  5 

DISTINCT  PERIODS  AND  PERCENTAGE  OF  REFLECTED  AND 
TRANSMIHEO  ENERGIES  PASSING  THROUGH  THE 
VERTICAL  BOUNDARY  (THE  CASE  WHERE  ALL 
three  zones  have  THE  BASIN  AND  RANGE 


PROPERTIES)  THE  LAST  TWO  COLUMNS 
are  the  result  OF  H.  MACK 

Period 

Sec 

Real 

Complex 

H.  Mack  (Complex) 

Reflected 

Energy 

Transmitted 

Energy 

Reflected 

Energy 

Transmitted 

Energy 

Period 

Sec 

Transmitted 

Energy 

10 

.0011 

99.9954 

.0073 

84.2120 

12 

.0054 

99.9945 

.0112 

88.2821 

13.8 

87.5 

15 

.0007 

99.9994 

.0017 

92.5926 

14.2 

94.0 

18 

.0010 

99.9989 

.0024 

95.3736 

15.1 

99.4 

20 

.0005 

99.9988 

.0018 

96.4662 

20.5 

99.8 

22 

.0004 

100.0000 

.0011 

97.0742 

22.3 

96.6 

25 

.0001 

99.9996 

.0002 

97.2620 

25.6 

93.6 

30 

.0001 

99.9996 

.0055 

96  6856 

28.4 

93.4 

35 

.0001 

99.9996 

.0165 

96.1306 

32.0 

94.0 

40 

.0001 

99.9997 

.0186 

95.8209 

39.4 

92.0 

45 

.0001 

99.9997 

.0155 

95.7497 

42.7 

94.7 

50 

.0001 

100.0000 

.0135 

95.7523 

46.5 

96.0 

60 

.0001 

99.9998 

.0122 

95.8963 

51.2 

95.3 

70 

.0001 

99.9996 

.0223 

96.0832 

80 


0001  100.0000 


.0001 


96.1831 
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TABLE  6 

DISTINCT  PERIODS  AND  PERCENTAGE  OF  REFLECTED  AND 
TRANSMIHED  ENERGIES  PASSING  THROUGH  THE 
VERTICAL  BOUNDARY  (THE  CASE  WHERE  THE 
VERTICAL  BOUNDARY  IS  AT  THE  RIGHT 
HAND  SIDE  OF  THE  ZONE  L 


Period 

Sec 

Complex 

Real 

Reflected 

Energy 

Transmitted 

Energy 

Ref 1 ected 
Energy 

Transmitted 

Energy 

10 

.1795 

1 

99.1060 

.2581 

99.7302 

12 

.0733 

99.3311 

.0551 

99.9465 

15 

.1209 

99.3781 

.1396 

99.8880 

18 

2.2780 

97.1127 

1.9161 

98.9499 

20 

.3447 

99.2575 

.3268 

99.6688 

22 

1.0353 

98.5935 

1.0766 

98.9225 

25 

.9698 

98.6958 

1.0128 

98.9857 

30 

.3175 

99.3742 

.3549 

99.6445 

35 

.4085 

99.2819 

.4607 

99.5385 

40 

.1525 

99.5212 

.1645 

99.8983 

45  i 

1 

.1209 

99.4230 

.1318 

99.8678 

50 

.0832 

99.4036 

.0892 

99.9100 

60 

.0555 

99.3044 

.0586 

99.9413 

70 

.1309 

99.1301 

.1081 

99.8914 

80 

.0454 

99.0328 

.0439 

99.9560 
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TABLE  6A 

DISTINCT  PERIODS  AND  PERCENTAGE  OF  REFLECTED  AND 
TRANSMIHED  ENERGIES  PASSING  THROUGH  THE 
VERTICAL  BOUNDARY  {THE  CASE  WHERE  THE 
VERTICAL  BOUNDARY  IS  AT  THE  LEFT 
HAND  SIDE  OF  THE  ZONE  R) 


Period 

Sec 

Real 

Complex 

Reflected 

Energy 

Transmitted 

Energy 

Reflected 

Energy 

Transmitted 

Energy 

10 

.0127 

99.9840 

1 

.0259 

84.2332 

12 

.0071 

99.9923 

.0123 

88.2958 

15 

.0048 

99.9943 

.0061 

92.5954 

18 

.0344 

99.9659 

.0342 

95.3359 

20 

.0145 

99.9848 

.0137 

96.4257 

22 

.0438 

99.9541 

.0406 

97.0042 

25 

.0497 

99.9490 

.0468 

97.1677 

30 

.0344 

99.9655 

.0417 

96.6803 

35 

.0653 

99.9338 

.0690 

96.1376 

40 

.0125 

99.9872 

.0240 

95.8200 

45 

.0109 

99.9888 

.0206 

95.6108 

50 

.0106 

99.9895 

.0162 

95.6320 

60 

.0118 

99.9879 

.0149 

95.8032 

70 

.1074 

99.8923 

.0241 

95.9536 

80 

.0283 

99.9718 

.0235 

96.1427 
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TABLE  6B 

DISTINCT  PERIODS  AND  PERCENTAGE  OF  REFLECTED  AND 
TRANSMITTED  ENERGIES  PASSING  THROUGH  THE 
VERTICAL  BOUNDARY  (THE  CASE  WHERE  THE 
VERTICAL  BOUNDARY  IS  AT  THE 
MIDDLE  OF  THE  ZONE  I) 


Period 

Sec 


Real 

Complex 

Reflected 

Energy 

Transmitted 

Energy 

Reflected 

Energy 

Transmi tted 
Energy 

.1410 

99.8621 

.1414 

91.7903 

.0710 

99.9305 

.0583 

93.9755 

.0537 

99.9441 

.0373 

96.2030 

.9384 

99.0587 

.8457 

96.3639 

.1468 

99.8515 

.1430 

98.0320 

.5724 

99.4253 

.5377 

97.9097 

.2310 

99.7686 

.2317 

98.3317 

.1871 

99.7836 

.1913 

98.1774 

.2673 

99.7323 

.2732 

97.8066 

.3002 

99.7002 

.3195 

97.5468 

.2449 

99.7547 

.2573 

97.4167 

.2027 

99.7973 

.2126 

97.4084 

.1677 

99.8321 

.1750 

97.6243 

.4799 

99.5201 

.4220 

97.1232 

.6412 

99.9581 

.0414 

97.5714 
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TABLE  6C 

DISTINCT  PERIODS  AND  PERCENTAGE  OF  REFLECTED  AND 
TRANSMITTED  ENERGIES  PASSING  THROUGH  THE 
VERTICAL  BOUNDARY  {THE  CASE  WHERE  THE 
POSITION  OF  THE  VERTICAL  BOUNDARY 
GRADUALLY  CHANGES) 


Complex 


Period 

Sec 

Reflected 

Energy 

Transmitted 

Energy 

Reflected 

Energy 

Transmitted 

Energy 

10 

.0177 

99.9815 

.0383 

87.122 

12 

.0119 

1 

99.9876 

.0165 

90.5295 

15 

.0092 

99.9906 

.0142 

94.0460 

18 

.0207 

99.9767 

.0266 

96.2619 

20 

.0073 

99.9912 

.0080 

97.1105 

22 

.0134 

99.9845 

.0084 

97.5348 

25 

.0271 

99.9707 

.0226 

97.5493 

30 

.0278 

99.9713 

.0114 

96.9659 

35 

.0577 

99.9379 

.0218 

96.3299 

40 

.0266 

99.9728 

1 

.0081 

95.9103 

45 

.0140 

99.9856  1 

.0065 

95.7905 

50 

.0111 

99.9886  i 

.0077 

95.7517 

60 

.0113 

99.9879 

.0107 

95.8399 

70 

.1399 

99.9598 

.0582. 

95.8825 

80 

.0209 

99.9790 

.0193 

96.0392 
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TABLE  7 

STARTING  COMPLEX  Vp  AND  Vg  VELOCITIES  OF  THE 
LAYERED  SYSTEMS,  THE  IMAGINARY  PART 
BEING  .n  OF  THE  REAL  PART 


Left  Lavers  j 

Rioht 

Layers 

i 

Vs 

Vo 

Vs 

REAL 

IMAG., 

(xlO-2) 

REAL 

IMAG„ 

(xl0"2) 

REAL 

IMAG, 

{xl0“2) 

REAL 

IMAG, 

(xlO-2) 

6.000 

0.  600 

3.540 

0.354 

5.800 

0.580 

3.470 

0.347 

6.000 

0.600 

3.560 

0.356 

6.700 

0.670 

3.740 

0.374 

6.700 

0.670 

3.600 

0.360 

6.700 

0.670 

3.740 

0.374 

6.695 

0.670 

3.595 

0.360 

6.700 

0.670 

3.740 

0.374 

6.650 

0.665 

3.550 

0.355 

6.700 

0.670 

3.740 

0.374 

7.900 

0.790 

4.430 

0.443 

7.000 

0.700 

3.820 

0.382 

7.900 

0.790 

4.425 

0.443 

8.050 

0.805 

4.620 

0.462 

7.899 

0.790 

4.420 

0.442 

8.050 

0.805 

4.620 

0.462 

7.880 

0.788 

4.390 

0.439 

8.050 

0.805 

4.620 

0.462 

7.860 

0.786 

4.355 

0.435 

8.100 

0.810 

4.540 

0.454 

7.820 

0.782 

4.320 

0.432 

8.400 

0.840 

4.660 

0.466 

7.678 

0.768 

4.182 

0.418 

8.500 

0.850 

4.500 

0.450 

7.730 

0.773 

4.220 

0.422 

8.500 

0.850 

4.500 

0.450 

7.800 

0.780 

4.278 

0.428 

8.500 

0.850 

4.500 

0.450 

7.890 

0.789 

4.430 

0.443 

8.500 

0.850 

4.500 

0.450 
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TABLE  8 


THE  CONSTANT  Qo,  AND  CORRESPONDING  Qa,  AND  THE 
CORRESPONDING  IMAGINARY  PARTS  OF  V.  AND 
OF  THE  LAYERED  REGIONS 


Left  Layers  (Og  *  400) 


IMAG  (Vq)  IMAG  (Vj) 
(xlO-2)  (xio-2y 


Rlqht  Layers  (Qg  = 


IMAG  (Vn) 
{xlO-2) 


IMAG  (Vc) 
(xlO-2) 
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TABLE  9 

THE  Qe  MODEL  OF  LEE  AND  SOLOMON  AND  CORRESPOND INfi 
AND  CORRESPONDING  IMAGINARY  PARTS  OF  V-  AND  V, 
OF  THE  LAYERED  REGIONS 


312 

672 

0.40 

6.00  i 

1000 

2095 

0.14 

0.17 

312 

664 

0.40 

6.00  * 

1000 

2400 

0.14 

0.18 

312 

810 

0.40 

6.00  1 

1000 

2400 

0.14 

0.18 

770 

2000 

0.20 

0.23  ! 

1000 

2400 

0.14 

0.18 

770 

2026 

0.20 

0.23 

1000 

2400 

C.14 

0.18 

770 

1836 

0.20 

0.23  i 

1000 

2518 

0.14 

0.19 

770 

1836 

0.24 

0.28  ' 

1000 

2277 

0.17 

0.23 

770 

1836 

0.24 

0.28 

1000 

2277 

0.17 

0.23 

770 

1825 

0.24 

0.28 

1000 

2277 

0.17 

0.25 

17 

40 

9.00 

13.00 

1000 

2436 

0.17 

0.23 

17 

40 

9. 00 

13.00 

1000 

2675 

0.16 

0.23 

17 

42 

9.00 

12.00 

1000 

2675 

0.16 

0.23 

17 

43 

9.00 

12.00 

1000 

2675 

0.16 

0.23 

17 

43 

9.00 

13.00 

44 

117 

3.60 

5.00 

34 

80 

5.00 

7.00 

44 

117 

3.60 

5.00 

Riaht  Layers 


IMAG  (Vd)  IMAG 


Left  Lavers 


IMAG  (Vp 
(xicr^') 
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TABLE  10 

FUNDAMENTAL,  AND  FIRST  AND  SECOND  HIGHER  MODE  EIGENVALUES, 
AND  CORRESPONDING  PHASE  VELOCITIES  AND  AHENUATIONS 
OF  THE  LEFT  LAYERS  (THE  CASE  WHERE  THE  IMAGINARY 
PART  OF  THE  INPUT  VELOCITY  IS  .U 
OF  ITS  REAL  PART) 


Wavenumber  (KM‘ 

•1) 

Phase  Velocity  (KM/SEC) 

Attenuation 

(xlO-Z) 

Funda- 

1st 

2nd 

Funda- 

1st 

2nd 

Funda- 

mental 

Hiqher 

Hiqher 

mental 

Hiqher 

Hiqher 

mental 

Mode 

Mode 

Mode 

Mode 

Mode 

Mode 

Mode 

0.2009 

0.1618 

0.1589 

3.475 

4.315 

4.394 

0.199 

0.1806 

0.1453 

0.1421 

3.478 

4.323 

4.423 

0.202 

0.1497 

0.1207 

0.1207 

3.497 

4.337 

4.481 

0.208 

0.1180 

0.0960 

— 

3.551 

4.363 

0.218 

0.0962 

— 

~ 

3.627 

— 

— 

0.226 

0.0853 

— 

— 

3.681 

— 

— 

0.228 

0.0765 

— 

— 

3.732 

— 

— 

0.227 

0.0662 

— 

— 

3.796 

— 

— 

0.224 

0.0542 

~ 

— 

3.866 

— 

— 

0.215 

0.0460 

~ 

— 

3.902 

— 

— 

0.209 

0.0400 

— 

— 

3.924 

— 

— 

0.207 

0.03  55 

~ 

— 

3.940 

— 

— 

0.206 

0.0318 

~ 

— 

3.955 

— 

~ 

0.208 

0.0262 

“ 

~ 

3.997 

— 

— 

0.216 

0. 0221 

“ 

— 

4.066 

— 

~ 

0.230 
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TABLE  n 

FUNDAMENTAL,  AND  FIRST  AND  SECOND  HIGHER  MODE  EIGENVALUES, 
AND  CORRESPONDING  PHASE  VELOCITIES  AND  ATTENUATIONS 
OF  THE  RIGHT  LAYERS  (THE  CASE  WHERE  THE  IMAGINARY 
PART  OF  THE  INPUT  VELOCITY  IS 
OF  ITS  REAL  PART) 
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TABLE  12 

FUNDAMENTAL,  AND  FIRST  AND  SECOND  HIGHER  MODE  EIGENVALUES, 
AND  CORRESPONDING  PHASE  VELOCITIES  AND 
ATTENUATIONS  OF  THE  LEFT  LAYERS 


(THE  CASE 

WHERE  Op 

IS  CONSTANT) 

Wavenumber  (KM 

-1) 

Phase  Velocity 

Attenuation 

(xlO-2) 

Funda¬ 

1st 

Funda¬ 

1st 

2nd 

mental 

Hiaher 

Hi  0  her 

mental 

Higher 

Higher 

Mode 

Mode 

Mode 

Mode 

Mode 

0.2009 

0.1618 

0.1589 

3.475 

4.316 

4.394 

0.2053 

0.1806 

0.1453 

0.1420 

3.478 

4.323 

4.423 

0.2053 

0.1497 

0.1207 

— 

3.497 

4.337 

— 

0.2065 

0.1179 

0.0960 

— 

3.551 

4.363 

— 

0.2092 

0.0962 

— 

— 

3.627 

— 

— 

0.2092 

0. 0853 

— 

— 

3.681 

— 

— 

0.2068 

0.0765 

— 

3.732 

— 

— 

0.2026 

0.0662 

— 

3.797 

— 

— 

0.1951 

0.0542 

— 

3.866 

— 

— 

0.1848 

0.0460 

“ 

— 

3.903 

— 

— 

0.1794 

0.0400 

— 

— 

3.924 

— 

— 

0.1772 

0.0354 

— 

3.939 

— 

— 

0.1768 

0.0318 

— 

3.955 

— 

0.1781 

0.0262 

~ 

3.997 

— 

— 

0.1836 

0.0221 

— 

__ 

4.067 

— 

~ 

0.1940 
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TABLE  13 

FUNDAMENTAL,  AND  FIRST  AND  SECOND  HIGHER  MODE  EIGENVALUES, 
AND  CORRESPONDING  PHASE  VELOCITIES  AND 
ATTENUATIONS  OF  THE  RIGHT  LAYERS 
(THE  CASE  WHERE  Qg  IS  CONSTANT) 
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TABLE  14 

FUNDAMENTAL,  AND  FIRST  AND  SECOND  HIGHER  MODE  EIGENVALUES, 
AND  CORRESPONDING  PHASE  VELOCITIES  AND  ATTENUATIONS 
OF  THE  LEFT  LAYERS  (THE  CASE  OF 
LEE  AND  SOLOMON'S  MODEL) 


Wavenumber  (KM' 

') 

Phase  Velocity 

Attenuation 

(xl0*2) 

Funda¬ 

mental 

Mode 

1st 

Higher 

Mode 

2nd 

Higher 

M^e 

Funda¬ 

mental 

Mode 

1st 

Higher 

Mode 

2nd 

Higher 

Mode 

0.2008 

0.1613 

0.1591 

3.476 

4.329 

4.389 

2.1463 

0.1806 

0.1449 

0.1421 

3.479 

4.335 

4.419 

2.0031 

0.1497 

0.1205 

— 

3.498 

4.345 

— 

1.7263 

0.1179 

0, 0959 

3.552 

4.368 

— 

1.3356 

0.0962 

— 

— 

3.628 

— 

— 

1.0034 

0. 0853 

— 

— 

3.682 

— 

0.8519 

0.07  65 

~ 

— 

3.733 

— 

— 

0.7852 

0. 0662 

— 

3.797 

— 

0.8549 

0.0541 

— 

— 

3.867 

— 

•• 

1.2567 

0.0460 

— 

— 

3.903 

— 

— 

1.7255 

0.0400 

— 

— 

3.925 

“ 

— 

2.1302 

0.0354 

— 

3.941 

— 

2.4612 

0.0318 

— 

— 

3.956 

— 

— 

2.7363 

0. 0262 

— 

” 

3.999 

— 

— 

3.1892 

0.0221 

— 

— 

4.068 

— 

3.6047 

83 


TABLE  15 

FUNDAMENTAL,  AND  FIRST  AND  SECOND  HIGHER  MODE  EIGENVALUES 
AND  CORRESPONDING  PHASE  VELOCITIES  AND  ATTENUATIONS 
OF  THE  RIGHT  LAYERS  (THE  CASE  OF 
LEE  AND  SOLOMON'S  MODEL) 


Wavenumber  (KM"^) 

Phase  Velocity 

Attenuation 

(xlO'2) 

Funda¬ 

mental 

Mode 

1st 

Hioher 

Mode 

2nd 

Hioher 

Mode 

Funda¬ 

mental 

Mode 

1st 

Higher 

Mode 

2nd 

Higher 

Mode 

Funda¬ 

mental 

Mode 

0.2002 

0.1543 

0.1534 

3.487 

4.524 

4.552 

0.0890 

0.1799 

0.1383 

0.1370 

3.493 

4.541 

4.586 

0.0891 

0.1490 

0.1149 

— 

3.514 

4.554 

— 

0.0901 

0.1173 

0.0915 

— 

3.569 

4.579 

— 

0.0928 

0.0958 

0.07  5  6 

— 

3.645 

4.615 

0.0958 

0. 084  9 

— 

— 

3.700 

— 

— 

0.0973 

0.0761 

— 

— 

3.755 

— 

— 

0.0986 

0.0656 

— 

— 

3.830 

— 

0.1012 

0.0533 

— 

3.928 

— 

— 

0.1176 

0.0449 

— 

— 

3.997 

— 

— 

0.1623 

0.0388 

— 

— 

4.045 

— 

— 

0.2346 

0.0342 

— 

— 

4.084 

— 

— 

0.3248 

0.0305 

— 

— 

4.120 

— 

— 

0.4256 

0.0249 

— 

— 

4.198 

— 

— 

0.6517 

0.0208 

— 

— 

4.308 

— 

0.9238 
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THE 


TABLE  16 

SIERRA  NEVADA  MODEL  OF  LYSMER  AND  DRAKE 
AND  CALCULATED  PHASE  VELOCITIES  USING 
LAYERED  SYSTEMS 


Thickness 


Compress ional 
Velocity 


Shear 
Vel ocity 


Density 


Period 


Phase  Velocity 


Lysmer  &  This 

Drake  Study 


6.0 
'.5 
5.0 
15.0 
20.0 


5.60 

6.00 

6.10 

6.90 

7.90 


3.10 

3.40 

3.40 

3.90 

4.39 


2.67 

2.70 

2.80 

2.80 

3.30 


6 

7 

8 
9 

10 


3.065  3.064 
3.105  3.104 
3.148  3.147 
3.191  3.191 
3.236  3.237 


85 


TABLE  17 

LEE  AND  SOLOMON'S  PHASE  VELOCITY 


Period 

Eastern  United  States 

Western  United  States 

Funda¬ 

mental 

Mode 

1st 

Higher 

Mode 

2nd 

Higher 

Mode 

Funda¬ 

mental 

Mode 

1st 

Higher 

Mode 

2nd 

Higher 

Mode 

10 

3.532 

4.468 

4.521 

3.379 

4.454 

4.657 

12 

3.559 

4.474 

4.543 

3.420 

4.511 

4.758 

15 

3.620 

4.483 

4.577 

3.489 

4.578 

4.980 

18 

3.697 

4.494 

4.618 

3.562 

4.658 

5.150 

20 

3.752 

4.502 

4.652 

3.611 

4.723 

5.232 

22 

3.805 

4.512 

4.691 

3.657 

4.799 

— 

25 

3.874 

4.528 

4.763 

3.720 

4.918 

— 

30 

3.958 

4.566 

4.926 

3.802 

5.090 

— 

35 

4.006 

4.616 

— 

3.860 

5.209 

— 

45 

4.047 

~ 

— 

3.941 

— 

— 

50 

4.054 

— 

— 

3.975 

— 

“ 

60 

4.061 

— 

— 

4.044 

— 

— 

70 

4.064 

— 

— 

4.116 

~ 

— 

80 

4.070 

— 

— 

4.190 

~ 

— 
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TABLE  18 

LEE  AND  SOLOMON'S  ATTENUATION  (x  10"^) 


Eastern  United 

States 

Western 

United 

States 

Period 

Funda¬ 

mental 

Mode 

1st 

Higher 

Mode 

2nd 

Higher 

Mode 

Funda¬ 

mental 

Mode 

1st 

Higher 

Mode 

2nd 

Higher 

Mode 

10 

.0903 

2.1459 

1.7652 

.2516  . 

.9800 

4.6070 

12 

.0933 

2.1146 

1.6871 

.2348 

.1986 

4.1221 

15 

.0958 

2.0682 

1.6130 

.2125 

2.9181 

3.9127 

18 

.0977 

2.0344 

1.6255 

.2017 

3.4195 

3.7966 

20 

.0982 

2.0239 

1.6693 

.2122 

3.6140 

3.6374 

22 

.0981 

2.0230 

1 .7286 

.2482 

3.7146 

— 

25 

.0983 

2.0344 

1.8300 

.3657 

3.7300 

30 

.1089 

2.0665 

2.0030 

.7190 

3.6524 

~ 

35 

.1446 

2.0970 

— 

1.1676 

3.6649 

~ 

45 

.2848 

2.1568 

2.0141 

— 

~ 

50 

.3737 

2.2036 

— 

2.3602 

— 

— 

60 

.5602 

2.3752 

— 

2.8826 

— 

— 

70 

.7409 

— 

— 

3.1902 

— 

80 

.9083 

— 

— 

3.3239 

— 

— 
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MAXWELL  ELEMENT 


n 


n 


Figure  2.  ELEMENTAL  MODEL  OF  VISCOELASTIC  SOLID 
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L  I  R 


Figure  4.  A  FINITE  ELEMENT  NET,  BOUNDED  BY  TWO  SEMI-INFINITE 
HORIZONTALLY  LAYERED  SYSTEM 


Figure  5.  COORDINATES  OF  AN  ELEMENT,  n  AND  C  ARE  THE  LOCAL 
COORDINATES  OF  THE  ELEMENT 
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PHASE  VELOCITY,  KM/SEC 
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FIGURE  7.  FUNDAMENTAL  AND  FIRST  AND  SECOND  HIGHER  MODE  RAYLEIGH 
WAVE  PHASE  VELOCITIES,  BASIN  AND  RANGE  PROVINCES 


PHASE  VELOCITY,  KM/SEC 
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5  10  20  30  40  50  60 

PERIOD,  SEC 

Figure  8.  FUNDAMENTAL  AND  FIRST  AND  SECOND  HIGHER  MODE  RAYLEIGH 
WAVE  PHASE  VELOCITIES,  CANADIAN  SHIELD 


PHASE  VELOCITY,  KM/SEC 
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Figure  9.  FUNDAMENTAL  RAYLEIGH  WAVE  PHASE  VELOCITIES  OF  ZONE  I 
(THE  CASE  WHERE  ALL  THREE  ZONES  ARE  ASSUMED  TO  BE  THE 
BASIN  AND  RANGE  PROVINCE.) 
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Figure  9A.  FUNDAMENTAL  RAYLEIGH  WAVE  PHASE  VELOCITIES  OF  ZONE  I 
(THE  CASE  WHERE  THE  VERTICAL  BOUNDARY  IS  AT  THE  RIGHT 
HAND  SIDE  OF  THE  ZONE  L). 
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16 


PERIOD,  SEC 

Figure  9B.  FUNDAMENTAL  RAYLEIGH  WAVE  PHASE  VELOCITIES  OF  ZONE  I 
(THE  CASE  WHERE  THE  VERTICAL  BOUNDARY  IS  AT  THE  LEFT 
HAND  SIDE  OF  THE  ZONE  R). 


PHASE  VELOCITY.  KM/SEC 
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PERIOD,  SEC 

Figure  9C.  FUNDAMENTAL  RAYLEIGH  WAVE  PHASE  VELOCITIES  OF  ZONE  I 
(THE  CASE  WHERE  THE  VERTICAL  BOUNDARY  IS  AT  THE  MIDDLE 
OF  ZONE  I). 
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Figure  90.  FUNDAMENTAL  RAYLEIGH  WAVE  PHASE  VELOCITIES  OF  ZONE  I 
(THE  CASE  WHERE  THE  POSITION  OF  VERTICAL  BOUNDARY 
GRADUALLY  CHANGES). 
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Figure  10.  RAYLEIGH  WAVE  ATTENUATION,  BASIN  AND  RANGE  PROVINCES. 

(The  case  where  the  Imaginary  part  of  the  input  velocity 
is  .U  of  its  real  part) 
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PERIOD,  SEC 

Figure  12.  RAYLEIGH  WAVE  ATTENUATION.  BASIN  AND  RANGE  PROVINCES 
(The  case  where  is  constant) 
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Figure  13.  RAYLEIGH  WAVE  AnENUATION,  CANADIAN  SHIELD 
(The  case  where  Qg  is  constant) 
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HORIZONTAL  VERTICAL 


Figure  16.  HORIZONTAL  AND  VERTICAL  AMPLITUDES  FOR  THE  FUNDAMENTAL 

RAYLEIGH  MODE  OF  SIERRA  NEVADA  FOR  A  PERIOD  OF  6  SECONDS. 
THE  CIRCLES  ARE  FROM  LYSMER  AND  DRAKE  (80). 
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Figure  17.  FUNDAMENTAL  AND  FIRST  AND  SECOND  HIGHER  MODE  PHASE 
VELOCITIES,  EASTERN  UNITED  STATES.  (Lee  and 
Solomon's  model) 


PHASE  VELOCITY,  KM/SEC 


Figure  18.  FUNDAMENTAL  AND  FIRST  AND  SECOND  HIGHER  MODE  PHASE 
VELOCITIES,  WESTERN  UNITED  STATES.  (Lee  and 
Solomon's  model) 
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19.  RAYLEIGH  WAVE  AHENUATION,  EASTERN  UNITED  STATES 
(Lee  and  Solomon's  model) 
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Figure  20.  RAYLEIGH  WAVE  AHENUATION,  WESTERN  UNITED  STATES 
(Lee  and  Solomon's  model) 


APPENDIX  C 


PROGRAM  RALEE 

This  program  is  a  revised  form  of  the  Program  Raleem,  written  by 
L.  A.  Drake  in  1972.  The  program  is  designed  to  analyze  the  motion  of 
harmonic  surface  waves  in  a  linearly  elastic  or  viscoelastic  layered 
medium.  The  layered  medium  is  subdivided  into  a  finite  element 
Region  I  and  one  or  two  semi  infinite  regions  L  and  R.  The  two  or  three 
regions  are  joined  together  at  the  finite  element  boundaries. 

The  program  consists  of  the  main  program  RALEE,  and  11  subroutines 
EIGVEC,  BOUNDA,  SECEVA,  SOLVE,  ORTHOG,  SETUP,  INVERTC,  MODAL,  STIFFER, 
WRDSC  AND  BANSOL. 

The  function  of  the  main  and  subroutines  are  as  follows. 

RALEE:  Reads  in  the  input  data,  assembles  the  stiffness  matrices 
of  the  finite  element  regions  and  calculates  amplitude  and  energies  at 
the  vertical  boundaries. 

EIGVEC:  Computes  the  entries  of  the  mass  and  elastic  or 
viscoelastic  matrices  of  the  layered  systems. 

BOUNDA:  Computes  boundary  matrices. 

SECEVA:  Constructs  the  eigenvalue  problem,  and  solves  it. 

SOLVE:  Solves  an  eigenvector  corresponding  to  t^e  eigenvalue, 
found  in  SECEVA. 

ORTHOG:  Finds  an  eigenvector  normal  to  the  one  that  has  already 
been  found. 

SETUP:  Decomposes  the  matrix  V,  used  in  the  boundary  matrix  of 
the  subroutine,  BOUNDA,  into  a  vector;  rewrites  this  vector  into  its 
real  and  imaginary  parts. 
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INVERTC:  Inverts  a  matrix. 

MODAL:  This  subroutine  Is  used  only  In  the  elastic  case  and 
selects  the  compatible  eigenvalues. 

SUFFER:  Computes  the  stiffness  matrix  of  an  element. 

WRDSC:  Writes  vector  B  and  the  matrix  A  of  equation  AX  «  B 
on  Tape  2. 

BANSOL:  Solves  equation  AX  «  B. 

INPUT  DATA 

Read  In: 

I.  Layered  System 

A.  Period  (F10.5)  (one  card). 

B.  Number  of  column  and  width  of  the  finite  element  region 
(110,  7F10.5)  (one  card). 

C.  Number  of  elements  having  the  same  depth,  young  modulus, 
rigidity  and  width  (6  110)  (one  card). 

D.  Number  of  natural  layers  and  depth  of  the  top  of  the  first 
layer  (110,  7F10.5)  (one  card). 

E.  Number, thickness,  compresslonal  and  shear  velocities  and 
density  of  the  sublayers  (110,  7F10.5)  (N  cards,  N  <  15). 

Step  D  and  E  may  be  repeated  once  If  there  are  two  layered 
systems. 

II.  Finite  Element  Region 

A.  Depth  to  corner  of  the  element  and  the  number  of  elements 

having  the  same  depth  (4(F10.5,  110))  (Number  of  cards  needed 
varies  In  this  and  the  next  two  cases.  For  example,  there 
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are  four  values  of  depth  on  each  card  and  each  value  multiplies 
by  the  number  of  the  element  having  the  same  depth.  Addition 
of  these  values  plus  the  values  of  the  following  cards  must 
add  up  to  be  the  number  of  elements,  N  ). 

B.  Moduli  of  the  elements  and  the  number  of  elements  having  the 
same  moduli  (2(2  FIO.S,  110)). 

C.  Width  of  the  elements  and  the  number  of  elements  having  the 
same  width  (4(F10.5,  110)). 

A  typical  Read  In  data,  described  above,  is  given  at  the  end  of 
the  program  RALEE. 
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0}»aH  ?a.P£  (TncjT, tapes, OJ7fUT,T4»El,rac£?) 


'Jl  lE 


,  .  H{10,  ?)  .Ml  J  1 


HH(ie) ,  . 
KKC 


(13,  15)  ,  Hi/ do,  Jo)  r  4*1=  (33)  , 
25)) .r (15,?) ,iS (15) 


’HA?i(3:) 


HE(.253)  ,HF(433) 

CC  IPLEx' ML  (i|  a  sTj';/ 15.15)  ;«<(“,  fl»  .0(15,’)  .HH  (15,2)  ,t  (is,  ?> 

1, VA(33,3-),  V3(3r,3:)  ,C*  (3C,3))  ,FO(?a,3i:),  iPtSC,?:) 

2. ZT(3l,3j}«M(3^«33) 


1!  1 


11 


PE^.T^:*) 


rCMnDN  nT.NF,  -n.  (’.MS.rM*-,:?!:,  IC.II.Al  (  EDI  ,  t«  (  'J  J  ,bZl3 
|?J)  ,:3(3j)  , ‘1  (3':)  ,'  2  (ZD  ,LD  (51)  ,'11/3(31  )  ,E  (?»■)  ,P1(3  J)  . 
3.J2(j.).\/l(3D),V2(3'),F0.D(3;),V(3’,3.),Ai(12C) 
C3'!;iO')/9ANtP.3/N'I.M>l,NL!f3LK,  °(  5  '  )  ,  S  ( 63 . 34)  ,  C  (  3  ti  .  3  J ) 

►ij  =  31 
^•^^=3  4 

‘'E4j(5,3)  T 
-O=*1AT(F10  ,5) 

EtAD  l.'I.WM 
'■CE'-AT  (110, 7?!?  .5) 

C"lE5A  =  S.2c3  53  2  717:i5a57' 

FET^^  111 

F jE-iar (♦jMa.  3ol'J‘*n;  width 

•••IN’  I.H.W^u^.OmET,! 

••Z-H.l 

=  EAD  2,iH,ia, IK.TL 
PD^'IAT  (4113  ) 

'"••INT  112 

-  =  j’.'IAT  (*:hdpI7D'|-,  aii,r  lAYEPE)  IT’JDTJ^E'^*) 

’EAT  ND.  OF  LAYERS  AND  DEFTm  dF  SUE'^ADP 
0  ID  J=l,2 
-EAD  1,N.,D 
f'iMT  12 

12  FD2'1AT(4X,6H.AYP'->‘,1DH  T  HI  C  K  )£  * '' ,  3  "  H  :0'1=. 

2  I 
K=1 

ro  11  1=1, Nt 

‘rAD  l,S*,THl,aLC,?)ET,fHO 

-•  T'JT  l,Ni,  TH”,  Alo,  3ET  .PMC 

Z=THr/N3 
PC  11  l=i,ni 
K  =  <»1 
HO (K. J)=D 
'-'=  JtZ 
H(<, J)=Z 

H‘t  {<.  J)  =  '*HD*(  ALP**2-2  .  *PPT**») 

P (<, J) =?H0*Fpt**2 
i:  f(<,j)»^HC 
!>;=< 

J'  =* 
f'Z»N»l 
Nr  ='1*N 
M1=N’-1 
N2=N3-2 
M7  =  ‘J3-3 


,£3(311,52(30  ,54(30  ,2: 
.  '=2(30,31(3 


VEL.  SMEAE  OE..  DEN'I 


. 

/  ? 
>  .t 
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r,:  P  !  =  ltN 

,rftr'''‘4r^:’TD'':D5SE9  O0:»Ti,  C)  1  ^“ES  !  IDM  AL  F..  DCI 5  I E  i  ,  'HE8? 
y/PLOCiriES’ANL)  3ENEITIES  C?  E'.EIENTb 
:E4D  93l«  nr(I»  .«(!)  tEMtlW  . 

9cii(Hr (I) .«(!) 

o:i  FO’.'IAT  (4  (F10.5»  110) ) 

K=0 

-■0  73Z  1  =  1. H 
Jt=K<(I) 

DO  7J2  J=l, JE 
K=K»1 

702  H;(<)=Hr(I) 

K=T 

PO  733  1=1. N 
:C  733  J=2.F‘ 

KsK*! 

733  )  i/(I.  J)*HS(K) 

••  E.SD  701.  (HEC.)  .KKCT)  .  j=l.I  A) 

r<isr  7.1,(HE(I).KK(I).:=1.I\) 

731  FC^'IAT  (2  (?-l3  .F.Il!:)  ) 

K=) 

: 0  734  i=l. lA 
J£=<<(I) 

'■C  7)4  J=1.JE 
<=<♦1 

714  H^{<)*S£(i) 

K=3 

’'0  735  1  =  1. N 
PO  715  J=1.K 
<=•<♦1 

73=  H. (I. J)*HF(K) 

•  FA*^  731.  (HF{^)  .K<(T)  .T  =  l,.Kt 

FilNT  731.  (H£(I>.KK(I)  .1  =  1.10 
Ks) 

•  j  716  1=1.1< 

JE»<<(I) 

'•0  736  J=1.JE 
<=<♦1 

73F  hF(K)*H£(I) 

F=3 

<■0  717  1  =  1. N 
'.•0  7  37  J*1.M 
<=<♦1 

737  D  (I.J)*HF(K) 

■iEAj  931.  <H!  (T)  ,«(T)  ,:«i,IU 
lul. (Hr (1) .KK(1I  .1  =  1. I.) 

K=P 

1 0  73F  1=1.1. 

JE=<<(1) 

PC  736  J=1.JE 
<=<♦1 

7ir  H:(<)*HTtI) 

K=1 

.0  739  1*1. N 
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CO.U''NS*) 


DP 


739  J=1.M 
K=<^1 

^39  :(I,J)=i3(K) 

17  1=1, Ni 
'0  17  Jsl.i 

K'.  (!  .  J)=*:  (I .  j)  *  (H-  (I»  J)**2-2«  J>  ••  2> 

17  r  (I, J) =<(I, J) *3 (i, J) 

hw(.M2  ,1>=HV(N,1)»H{N,:) 
k-VCN:  ,1Z  )=HV(NZ  ,1) 

:o  IB  j=2,M 
HVINZ  ,J)=HV(NZ  ,1) 

1*  HH(J)  =  (J-1) 

HH(1)=3. 

HH(MZ) sWB 
P'TNT  2j 

?r  rcR14T(*CX  CDO^Cl NtTEf  OP 
PnINT  22, (HH( J) ,J=1,M2) 

2?  (ix,21-P.  2) 

Jr-.I-JT  233 

23r  FO^BAT (•uSEPTH  TO  OO’NE^i 
:0  2i»  1  =  1. KZ 

21*  ^^Tn:  22.  (HV(I,  J),  J=l,PZ) 

-RIN’-  229 

Z?*"  F07*«4T  (♦  jLAKOft  CP  ELE^iENTE*) 

’■J  2=>  1  =  1. N' 

2F  PRI'JP  22.  J)  .J  =  1.’M 

FRIN7  213 

’f  PO-^IAT  (♦CiHEiR  “IJOMLUS  OP  ELiiEN^s*) 

■.0  21  1  =  1.  N 

21  22.  (3  (I.  J)  .  J=l.f^) 

F-tlMT  22  3 

’7''  P0".'14T  (♦3DENEIT  Y  OF  ELF*1EN'r*» 

•0  25  1=1, N 

2’  FV.IN'  22,  (’ (I ,  J) .  J=l.  M) 

C‘-:*:=0ME3fi*0*lE06 

C  3FNEP4TE  £QUI\/A.E*^T  Br.'JNOARv  vAfilX 
[0  57  1=1. N 
t-A(!)  =H(:.2) 
hL  A(I)  =H'1(1 ,2) 

GA(I)=^(I,2) 

57  •-£(I)=S(I,2> 

:A..  EliVECC-iA.^^-A.GA.f  A,N) 

^RINT  It  01 

1321  format  (•CElfcENrfAuUEP  OF  Tr.A  n;'1’'7TE0 
eRI'J7  59, £ 

513  PC»MAT  (IX,  1CE13 .6) 

FF=E(IR)*(3.,1.) 

CC  ?51  1=1, 

JFfPEA'-(E(I)»  .-P.''.>  GC  TO  3il 
V\/=01PGA/5:£A.  (E  (’  ) ) 

C=A11AG  (£(!)»  /REAL  (Ed)) 

3F(Q.EQ.3.)  30  7j  3*^? 

0=0*2. 
frT!^7  55,  g 
XRX  PRImT  59,Vtf 
551  CO!)I:n'Je. 

:<  =  IR 


-iOOET*) 
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CD  61  J=1,NP 
e.  f Al J»*V< J.Ift) 

6!  ‘"TCJ)=£{J) 

F^pU  59.  (tf  <I,I'),I  =  l,rP) 

^AL.  90UNDA  (HA,  -■'A.HLA.GA.L" 

P1  =  J  • 

P2=0. 

'••'J  723  1=1,  N» 

P?(I)=0. 

PC(I)=1. 

'0  745  J=1,N3 
F3(I)=P8(1) ♦W(I, J)*PA(J) 

>  FC(I )=p: (I) ♦AIhAG (W (T, j) ) ♦3A( J) 
P1  =  -H'DCNJD(PA{T)  »*P5(:) 

!  P2=P2^C3NJ3(PA(I) )*3C(:) 

FlsOMESA/Z.’AlHAGtPl) 

Q=?2/Pl 


pcTNT  59,0 
CO  211  1  =  1, 

:0  231  J=1,N3 
ZT(I, j)=w(I,J) 

?D1  V5(I,  J)  =  i/(I,J) 

:  35N5?A’’£  E3J1  3DJNr;F,Y  MA‘-'».TX 

CO  56  1=1, N 
HA(r)=H(i,i) 
h.A(I) =HH(I,1» 

CA(I)=r(I,l) 

56  -AtTjs^d,!) 

CflL_  FI3v(ED(HA,H.a,3A,rfl,N) 

PtTnT  1003 

ID*??  KC'>'<4T  (•■jEIGENVAlJES  pr  iNDDENl  130 
59. E 

FCs£(I3) *(:.,!.) 

C'O  362  1  =  1, NP 

IF(^EAu(E{I)> .LE.O.)  rc  TO  3>2 
VV»D1EGt/REA. (P ( r ) ) 

C=AI  lAGCE (i ) ) /PPi. (E (T ) ) 

IF(Q.EQ.j.»  30  TO  359 
0=3*2. 
fp.TMT  59.  Q 
(6«»  5S.Vt( 

35?  OC'ITIN'JE 

14  =  T1 

PC  222  1=1, NP 
?2?  PP(I) *rf(! ,lw) 

Ps^NT  59,  (3K(I)  ,:«l,rjP) 

:  ptn-:'  3T3P' aceheit  »irpix  pop,  9F-.pr'P3 


ro  63  !«2,NP,2  --I  -j 

CO  53  J*1,NF 
vd,  J)  =•;(!,  J) 

L  p*0 

CAL-  90J;<LA  (Ha,PA,HLa,GA,L?,PD,W,=D, 

F1  =  J  . 

P2=3. 

CO  735  1=1, KP 


lODEi 


SytafaS^^&S 


JX3  PAGt  Li>  S*Ui 

K«  c^i  n>  u,a  _ ^ 
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5t' 

2  ”'2 
1! 


51 

“al:  j 

?7 


la- 

c  ■<<? 

3S 


r=(T»=T. 

tC(!)=!}. 

755  j=i,N3 

Pb(!»  =9^(1) ♦rt (I , J»  *9? ( J) 

FC(H»o:(!)  J) )  *?sc  J) 

Fls»i»C0NJ3  <39(T)I*  =  9<T) 

'■2=?2*35:JJS(9’.(I)  )*»C(I) 

»ls3iE5A/2.*AI'1ftG(Pl) 

Qsoz/pi 

c<IMT  53tQ 
•^C  51  1*1, NP 
[•0  51  JsI.NP 
W{J,I)  =  -I(I,J) 

3(J.!)*C(I. J) 
r-O  202  1*1, MS 
:.'0  202  J=1,NS 
V4(I,J)*V(I,  J) 
rc  11  1*1, NP 
GUI)  =  9P(I)  •inC  . 
no  91‘  1*2, NF, 2 
K=I-1 

'.C  91  J=2,NP,2 
L*  J-1 

V(T,J)*  -C(l.J) 

7(K,  J)  s  0(K,J) 

V(K,L)*  -C(K,L) 

V(I,.)*  C(I,.) 

.A’;  sTIs-N^S?  am:  CO'JilPTE'JT  ifl?:  1AT5..CE5  0-  SLOPING 
MHsMN^NN 
:C  27  1*1, NH 
"(I) *c. 

no  27  J=1,M) 

A(I, J)*C. 

*£<<fsD  2 
NU'19';.K  =  U 

I  ooj  n*3 

1W  =  G 
no  17C 

no  170  1=1, N 

Zrtu-Hti)  95,  •5,19r' 

IV*I7-N)4 

fE  A  Block  o-  EO'jA’roKT 

CAL'.  WxOSC  (LCCJNi  ,N9) 

10*17-1 
UUsfrf^l 
L7»LJ»1 
JU*IO».7= 

JV*JlJ*l 
KU»JV^1 
KV=KJ»l 
Y(l,l)artH(J) 

Y (2, 1) *HH ( J*l) 

Y  (3.  1)  *flM  (  J^l  ) 

Y(A,l)sHH(J) 

Y  (1,2)*H/  (I,  J) 

Y«2,2)*'17(I,J*1) 


(6U)  OM  1 ftaE  2 


7 

I) 

£ 

^  a 

->  -T» 

as 

4  9 
^  % 

••  >Nf 

•fir 


118 


Y<3,2)=M\/(Itli  J^l) 

Y(^,2) =HY(I^1,J) 

CfiL.  STlrF-.(Y,Hl.  C,  J)  ,f  (I,  J),  ■;  (I,  J)  .3M?,aK) 
:^(:.E'J.N)3C  TD  lOD 
i  (IJ#Nr'»3)  =A(IJi‘l»*3»  ♦/iK(lf5» 
a  (IU,N3*4)sA(I  J.NOtt*)  ♦*<(!,  6» 
fi(IJ.I)sA(rUt3) ♦&<<!, 7) 
a(lj,4)=a(ij,<*»  ♦ftKd.F) 
fl (I  si (I V.NS»2) ♦fiK  (2.5) 

a(IV.N?*3)sfi(lV.M='f3)«ftK(?,i) 
a(Itf.2l =A(IV.2» ♦4<(2f 7) 

4(I\/.3)«A(IV.3)  ♦fi<(2,e) 

A(  JJ.7) sA( JU.  3)  ♦AKiB, 5) 
a  ( JJ,4) =A  C JJ«  4)  tAK (3t  &) 

A ( JV.2) sA ( JV,  2) ♦AK(4,f ) 
a  <  j'/.3)  *;  ( JV.  3)  ♦a<(4, 6  ) 

A(<u.i»=a(<u.i)  ♦aKts.E) 

A(<U,2)sA(<U,2> ♦aK(^,6) 

A  (l<V«  1)  sA  (KV.  1)  *AK(n,  6) 
a(.J,l)=A(.U,l) ♦AK(^.7) 

«{LU.2)=A(.U, 2) ♦AK(7, 6) 
fi(LJ.N  =  -l)  =  a(L'J.NP-l )  ♦AK{3,  71 
A(!.U«NP)=AlLJ»'f  =  )*A<(4,7) 

A  (;.U.N?*1)=A(LJ.  •)»♦!»  ♦rKt5.7l 
A{LU.N“»2)=A(LU.'JP»?)  ♦AK<6.  7» 
Ali.V.l)=A(LV.l)  ♦A<(«l,e) 

A  (!.V,NO-2)  =  A  (LY  ,?;=-2)  ♦AKtS,  6» 
A(LY.S»-1)=A(LY,NB-1) ♦AK(4,  e> 
«(-V.NP|=A(LY.Ne)fa<<{D.E) 
A(*V.NP^l)*A<LV,:i®tl)^£(<(6.6> 

I"*"  A(TJ,1)  =A  (lU,  1)  ♦AK(1, 1) 

4(IJ,2)sA(lU.2» ♦a<(l,2) 
A(IJ,NP+1)sA(IJ,n3»1) fAKCl. 3) 
t(TJ.NP>2)=A(:j,,-J^>»2>  ♦fKd,  4t 
A  (lY.l)  =a  (I  V,  1)  ♦A'<(2,?) 


14? 


I5r 


71 

14? 


A(tV,NP)=A(IY,NP) ♦A<(2.3) 
MIV.MS*1)=C(TV,*IP»1)  ♦fK(2,4) 

A( JJ.1)=A(JJ,1) «a<(3.3) 
A(JU,2)sa(JU,2)  ♦AKCS.A) 
a(JY.l):A  ( JV.  1)  «^A<(4.4) 

:2=i*i 

.1)3C'  TO  141 
mf»nip-ii>i 
f::*NP-I2»l 
*■•2  14C  Jl^l.SF 
N&»ii^ji-i 

Adut Ji) >A(iJ« Ji) »:  (ii*NG) 

:o  15:  ji=i.si 

NC-»I2^  Jl-l 

A(IV.J1»»A(IY,J1> *C  (I2tNG» 

•?(IJ»=J. 

T;  ( !  V  )  s  0  . 

"O  71  K»1»NP 

=  (TJ)=5(IU)  ♦(r  (i  'J,K)-V(lLi.  KM  •«-(<) 
-(IV)  «3(IV)  ♦(C  CV.K)  -VdV.Kl  )  'c^CK) 
TP(J.Lr.'l)  3C  TD  I7r 


119 


145  Jlrl.MF 
No=Il» Jl-l 

a(jj,  J1)=A( JJ, Jl) -CT (I1,NG) 

•  0  155  JI  =  1,'1T 
NS*I2»Jl-l 

e  ( jt/.  jn  =  A<  jy,  ji)  -cr  (r?,NG) 

co'in^'jE 

'“ALL  W503C(L:31IMT,M^) 
i"  (-^OJNT-'J'!)  19  5.  ?1 5, 3 3: 

315  LG3Li‘n  =  LC0yNT*l 

~\rfjr  146, NUMB.  < 

FO’^'IAT  (*iNO.  OF  SLOOK^  OF  EQJ1’ID^5  =*.IU) 

'-COslCOJUT -1 
147, .CO 

147  FC1HAT(*  MO.  0"  tOja^yrM?  =*,I5) 

'TLVE  THE  EQUanrsr  WtTTTtH  O'J  '1-E  2 
cal.  5AN5C. 

P^.TNr  631 

FC»''AT  (•CMAXIHJh  (COH'I  cX)  DZi^L.  OF  CO''''JE'  »rTM*3*) 

no  63C  K=1,M^ 

FKOsK 

=  =<I  JT  6:2,  (K1  (I )  .  3fT)  ,  J=1  ,'JPl 
'■C.'lAT  (4(*5,2r  13.6)  ) 
f  AH’MrjTEi  AHO  SMa^F  cfi 


145 


15c 
IT*: 
1  •  s 


33: 

146 


631 

63' 
63  2 


)F  2E.-Li:-  EO  ►'D0E3 


CO-*=JTf  AH’lirjrEt  AHC  3Ma3E  .EAOE  Or 
CO  41  1=1, NF 
41  =(I)  =B(l) -05(1) 

F'nlfJT  2C,16 

2ni6  FO;  iAT  (•0\E"lECT -r  ‘1C'’E  n *  >L ACEHcN^ •  ) 

F^IHT  35,  (Bd), 1  =  1, HO) 

00  639  1=1. N3 
At(I)=5. 

CO  639  <=1,N» 

60'=  AA(I)=4A(I)4  tfa(I,K)*P(K) 

00  533  1=1, N3 
AMO(I)=:aB5(AA(I)) 

63?  PnA^EdlsATAN^lAiHAGlAadll.irfiLCAAd)  >) 
F”,rHT  6w4 

634  F02HAT  (♦!  AHP.TT'j- -r  ^  =  F.F:T"0  HOD"**) 
P'lVT  59, AMF 
F'.INf  636 

606^POPHAU''aPHA;E  LEAn'  or  p.EFL-IC'EO  HOOEo  (A 

prfHr'^59, phase 

HJ  =ieOC3 .7=>Ea'.  (E  (TH)  ) 

,  :0  514  1=1, n 

61'4  P(t,l)«4HF  d)  ••2/^EA!.(f  (D) 

F^TNT  623, HJ 

623  FO^HAT (•ilNCIOEH-  EHE9CY  =*,:i3.5» 

621 

621  Fo»-AT(*  9EF,EC’’£3  ENPrCY*) 

59,  (»<!,!>  ,1=1, Tf  ) 
r*Jl=HM»l 
J=J 

'U  616  IsHMl.NP 


L.H.  EO  OF 


.OPING 


120 


!1  ^ 
’017 

61  r 

61  •? 
63  5 

6?  7 


511 

615 

62? 


jsjtl 

FO*1fiT(*o  FIGH"  BO'JNne^Y  SP)E  riS’LM  F-?1i  V  •> 
FvINT  59. (Ss(I) .1=1,NP) 

'.0  612  1  =  1.  M’ 
i6(T) =3. 

LO  612  K=1.N’ 

®A(I)=AA(I)»  V3  (I.K)  *3  (K^KM) 

60  613  1=1. fP 
AMO(I)*CAE'(AA(T)l 

FHASE(I)=ATAn2( AI-AG(AA(I)) . iE AL ( A  A (I )  ) ) 

F'lNT  6d5  • 

rO’MAT  (♦0A'1P_ITU:Er  OP  Tp-AMSII ’■"53  MDDE>*> 
PP.INT  59.  AMP 
P».IMT  6£(7 

FO’VAT  (♦CPHA3E  LEAD'  OF  TPAN;mIT’’E3  IOp-Es  (AT 
2'-T^'lJCTU^Ei  *) 

=;INT  59, PHASE 

PA  (!)■=  (iii!  (ET  ( I) )  *Wi»PHA?P  C  )  )  /OHESA 
'/WsHi/PA  (I ) 

P  INT  59,y/V 

Z'j  615  1  =  1,15  .  , 

F  (!,2)  =A1P(1)  ♦•2/’EAl.  (FT  (I)  ) 
ff.INT  522 

FO’HAT  (*01  ^AN?httte:«  EKEPr-y*) 

FPIMT  59.  (P(I,2)  ,I=1.I'') 

P1=J . 

F2=3. 

P3=D. 

F4  =  3. 

F5=j  . 

P6=0. 

f.C  47  1*1, MP 
FA(I)=3. 

P5(T)=0. 

F5(I)=3. 
no  49  J*1.NF 

F;(I>=PA(I)  ♦W(!.J)*9?(J) 

P-s(I)=os(I)  ♦*((I,J»*’(  j) 

PC(I>=PC(H  ♦ZT(I,  J)*3:  (j)  . 

Pl  =  Pl+ REAL (6» (!)) ’AT  hag (PA<!»  » 

F2s?2>AIHAG (3R( !) » ♦peal (PA( T) ) 

F3=P3»R£AL  (E(I)  )*ArHAG(PB(:)» 

F4*»4^Al!-tA3  (3  (I  ) )  *RFAL  (P=  (I  )• 

''5=s54R£AL  (B>  (I)  )  •ATMAf.(OC  (  7l  > 

F6=?54AIMA5(35(T)  )*9rAL  (pcdi  ) 

EI»(Pl-»2) •.5*0HEGA 
£•’.=  (’3-54)  •,5*0HEGA 
ETT=(P5-P6» •.5*0H5GA 
EI*A3S  (EI) 

6>'.*A3S  (£R> 

FfsABSLEM  ) 

P’THT  3333 

3?T»  =ORiAT  (  •0If-3IDENT  ENERGY*) 

FrlMT  59. FI 


7.H,  EMD  0’ 


49 


47 


?.D»: YG 


121 


44  44 


55  55 


4«  e; 
3''7 

20  5(? 


4444 

F0-''*1AT(*0!’.E*LECTEn  ENERGY*) 

59,  E? 

F-slNT  5555 

F0?'1AT(*Cre.ftN3MlT’'ED  FNERGYM 
F^I^r  59,ETT 
Pl=£I*ETT 
F2  =  0. 

no  307  1=1, N» 

F:(!»=3, 

rz-  405  j=i,N» 

FC(:)*p:(I)  ♦(50  (i,  J)  -Ft'd,  J)»  ♦6F(  J) 

='2  =  S2*C0NJ3  {3R(  !)  )  *30  (f ) 

FnINT  205C 

FO?'1AT(*OCH4N3E  IN  HAVFHL'HEE^  EETi^EEM 
0K=?2*0'iE&6/(4*Pl) 

FRINr  59, CK 
F^lsETT^El 
EI2»ETT-£I 
:K=ET2/£ri*E(:^) 

P<TNT'59,rK 

f  Ml 


'^3  FOL'MDARIES*) 


122 


icri 


SUl%aUTl;j£  EIG>/EC(H,HL.G,P.M) 

C  H(15)  (15) 

COii^'.EX  c,EO.D,?l»^2,'Jl,U2,VU  V2.I/.A9 

‘c6r.:l5  Si:S3:5^:ir;3rJ5;l!(:ir-’''’''"='-‘"='-'='- 

Ca'iPlEX  HL(15),G(15) 

COinf^  '(0»NF,Nl,N2,N3,r*lSt!f,IC.I!.Al(  33)  ,C3(3j), 
23t,)«^3(30)  t5:i{2C).E2(3r),LD($J).M((3(30)tE(3C),nl 
3.y2n0)  ti/1  (33) ,  V2(3  3)  .TOLD  (31)  .V(33,30).A5(120) 

9  I=ltN® 

f  EOL?(r)=3, 

►'40=3. 

A  1/0=3. 

BHOsO. 

CH0=3. 

rvo=3. 

:  =  3 

DC  22  J=1,N 
H3=H ( J)/3. 

Ai.=H,.  (  J)  ♦&(  J)  ♦G  (  J) 

CMs2  ( J)  •H3*0’1S 
AH=A.*43 
Arf=G( J)*H3 
CH*G(J)/H(J)-C4 
CV=Al/H(J)-C4 
?«*(5( J)-H. (J) ) *.5 
:’V*(3(  J)  ^HL  (J)  )  *.5 
I  =  I^1 

£l (I ) =AM*AH0 
•*•3(1)  =A-(*.5 
E2(I)=3f1-E40 
cAd)  =9(/ 

rid)  =CH»CHO 
C’C)  s-CH-l  .5*C‘' 

Isl^l 

A1  (I)  >Ai/«'Ai/0 
A3(I)=Aw*.5 
‘■Ed)  =-5V 
6<.(I)  =  0. 

Cld)=C»/>CVO 

C3(I)»-C»/-1.5*C^' 

AHO»AH 
AVO=AV 
=iH0=9H 
THOsOH 
2?  CVO=CV 

' ''l.^|l^1^|,gDN-Ll4EAi  EIGFr.i/ALJF  ?'03.E4,  (E**2*A  ♦  E 
K=AI4A§(Hl) 

!F(K.EO.O.O)  CALI.  MOOAl 
EW? 


32(30  »3i*(33)  ,:i 
33),;;2(30t  Jl(3:. 


•3  ♦  C)*»/  =  C 


-0 


■.J  2u  J=1,NP 


t'U  J=1,.NP 

T"(.;<.-..T.i»r,o  T.1 
0=  d(J)*(0.,l.) 
3JT0  13 

12  :;  =  -£  (  J)*  (J  .  ,1.  ) 

13  JT  23  1=1, .4P 

20  :(i,j)=o(i,j) *0 

s>.7u? ( 7,no) 

=  I,NO 


12 


Jj  4.0  l! 
j  J  3  0  J- 
71 ( J)  =.. 
jJ  30  <=1,  4P 

71 { J>  =  V1 ( J) ♦C ( I, <) ‘VCK, J) 
J=  I, -JP 
-(I , J)=V1 ( J) 

■/  I  =  t  . 

YI  =  -. 


«=  J 

J3  5-  T  =  2,:42,2 
<=I-1 

•I  =  M  ♦  1 

X J= ."•HL (  1) 

Y  J=.5*b(‘1) 

-  (<,  I)=0  tXI-XJ 

J='<«’3 

J:I>1 

-(I,J)=C(I,J)>YJ 
XI  =  X  J 

so  Yisvj 
.i  =  .U  1 

-  (  ^^l  ,UP)  =v.  (N1  ♦  VI-.  '.'.•ml  (M) 

Is 


470=0, 

371=0, 
0  H  .3  =  0  , 


33» 

4*(30) 

G)  ,"1(3 


•)=  r.,  -  - 

»i  =  o. 
vi  =  -. 

:3311  I=£,rj?,? 

<=I-1 

IS 

♦J(H) 

:  1=^(*1»*H3*0riS 

-r(=4L*M» 

'f  i=-CMHHI-G<MI)*,9 
»  (HLCJ1)  ♦SC‘11  »*.? 

-i-CKfOs 

“<  U<,I>S  JHJ-3H 
*<{<» 

.;.*(;<, 5 
-•I  =4H*.  ■ 

(  I»1,K) =AH*.5 
w:(K,I*l)s-CH-l,5*C« 

H(<,K  +  3)='»H 
W1(<,<^3)S  3H 
W(  I,  I)  =A»/tAVO 
s-IA  ( I  ,I»  =Mi/*fi»/0 

n'-,  ( i  »r)  5wV>Ctfo 

tn(i.i»i)  =  Yj 
W(I,l*i)a3V 
(  I  t  :♦!)  s 
.M  I 

I«-2»  =  av,5 
-<(I»2,1)  =-\/*,5 
!U(I»2,I)=AV*.^ 

.I»2)=-CV-1.5*Cl1 

H  .■)  s  A  H 
4./0S4V 

':H2sr.M 
j  /nsctf 
Yl-syj 
XlxXJ 

;t=-(L('!)t3('i)4’bC3) 

2'1  =  ><(;l)XH3*Or3 
A  3*  AL*H3 


125 


i  I-  G 

3  7:  .  =  (Ml 

(•■1) 

3H:G  (*<)/H(M) 

XJs- (He (  1) -G (H) ) •  .5 
Y  Js  (HL(  IJ+GC?-*  )*." 
i!(i4l  ,:<!)  =4H  +  A‘^0 
WA  (  Ml  ,:il)  =4H-f  AHT 
^G(U1,M1)s:h»CH3 
^(MltMP):  (Xl-XJ) 

-((Nl  ,NP)  =  3H0-.5*HL<*I) 
.J3(Ml,M?)s3HO-.=i*HL(H) 

=AV  +  A70 

(  Nl-  ,ij?)  =47'*-A  V3 
A  Z  ( iMPtfi?)  iC7  +  GV3 
jO  413  1=1, MP 
Jj  ••liv  J=I,  ?iP' 

.(  3  (  J  ,1)  =W3(I,J) 

W.](J,1)=MD(I,J) 

P3IMT  =9,PC,FF 
^.^.•■.^'AT(lX,lOl  13.6) 

:f(l.«.i3.1)  &3  TO 
j3  til  J=1,MP 
nj  til  1=1, UP 
7(i,j)=w(i,j)*f; 

.-Ij  (  I  ,  J)  =  K3  (I  ,  J)  ♦FC 

-<=FG-Fe 

■!J  412  J  =  1,N? 

IJ  412  1=1, M3 

!j A ( i ,  J)  =  (1 ,  J )  •  cn+w3  ( * ,  J )  +wn  (  I ,  J) 

03  TO  5036 
GJIiTiNUt: 

■;<.  rp^'^F 
)J  -Cs  J=1,NP 
:  )  t09  1=1 ,NP 
.1(I,J)=K{I,J)»(-FF) 

-iJ{I,J)  =  ;0(I,J)*  (-FF) 
n(I,J)  =  KA(I, 

_J  'i-.  3“  I  =  l,Nf' 
jJ  3C33  J=1,M3 

44  (I,  J)  :  W4  (T  ,  J)  Y'WOCI,  J)  Y-mG  (  I,  J) 
i’p 
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CJ3^3lJTI.4E  StCEVfi 
CO*<P-EX  W(3C.3S) 

COMMON  ND,NP«N1  ,N2,‘J3,rMS,IP.IC.!I,4l!(  30)  ,  fi3  ( 3:  ) ,  B2  ( 3  u)  .  3U  ( 3  } »  ,  21 1 
230),:3(30)til(3C)tS2(3C),LO(J:)  .-rfSCSO  ),E(3:)  ,^1(33),P2  (3a)f  JUBO; 
3.U2(3C)  .</l(33)  «  V2(30)  •F0Un(3))  «  V(30.39)  te3(12C) 

DATA  EPS/I. CE-7/,EP?;2/l.CE-lJ/ 

NNl=N9-l 
NN2=MP*2 
00  73C  M  =  liK'Ml 
MV3(N)=)4*3*N0 
I  =  N»1 

700  ^<V3(I)=1♦2•N^ 

NV3(NN1) sNNI^NO 
XN  =  0, 

CO  15  J»1,MP 
XsCK  J)/Al(  J) 

15  XN*X)UA^r  (X) 

XN=2QRT(XN/N») 

r=CMPLX(XN,XM»XN) 

rc  132  N=1,N3 

LO(N>>0 
tr).2(N)=0. 
ri(N)s(u.«c.) 
t^ECN)  =  (0. 

•■KN)  =1. 

E2(N) *C 
VI <N) =S1 (N) 

10?  V2(M)»32<N) 
fVsXN/NP 
Xl«AlHAb(Al) 

IMXl.NE.C.C)  GO  TO  8Cr 
'.0  106  N=1.N» 

‘2(N)sl, 

108  V2(N)s3?(N) 

EV*!. 

on  900  N»1»NM1,2 

“V8(M>»Nt3*ND 

I*M^l 

^O*"  KV9(I)xI»2*nO 
*'V8(MN1)*NN1*N3 
80?  K*0 


HC*3 

IM«0 

rVjNT  2i3C 
P^INT  2037 


l*0C  K«<»1 


110 

8^ 

88 


IFfXl.NE.C.b)  GO  ro  88 
JF  (cDLOiO  .EQ.O  .:)  33  ' 
tV«E3L3(<)*(l.t0.30?l) 
GO  T3  87 
XsCA5S(£V)/2. 
ev*;mf.x (X  ,x) 
rD«i . 

DC  112  Nsl.NP 
U1{M)*31(N)*V1(N) 


3  lie 


127 


11? 


11 : 


'->5 


■'O 


75 


85 


8® 


U2(N)  =M{N) 
uC  113  Ssl,N2 

U1(N)bJ1(N) ♦:3(N)*V1(L} 

L'l  (  -  )  ='Jl  1 1 )  *33  (  N)  *71  (  N) 

U? (N) ®U2 (N) ♦ft  3 ( N ) *72  (  L  ) 

ui(i.)*'J2(L)*63(N)*72(N) 

iFtXl.cQ.C .C»  SO  TO  89 

IFd'l.iQ.l)  so  ^0  65 

IM*5 

00=1. 

X=Cfi95(E7) ^2. 

EVsCM^LX (X»  X^X) 

SO  TO  99  . 

£7=EV*(l.t .0001) 

!H=2 

CPsl. 

T  KsJ 

ro  230  IT=1,20 
FVS=EW**2 

72  ( N)  »  (U iH )- U2  <Ti *00 

NQ*M  ♦NO 
rif<sNQ*i4U 
HS  =  N’.*N0 


AS(NI 
ft&TNQ) 

AB(N>) 

12C  AQtNS) 

CEs(0..0.) 
CAls.  >Ol7E 
ro  llu  .N=1.N® 


=E75* A1 (Nl ♦Cl (T ) 
=  E7*32(*I) 

=E7E*A3 (N) ♦CS (t) 
sE7*94(N) 


^0 


;  '3 


:35 


IP 


13: 


172 


13: 


71{N)*(72(N»-71  (M)*C'')/E7 
CE  =  :Uil(N)*71(M)  ♦U2(N>*72(M 
r.E  =  CE*C0 
CO  132  NslfN® 

U1  (N»  sOH  N)  *71  <  N) 

02 (N)*A1 (N) •72(N) 

00  133  N=1.N2 
L=N^2 

01(N)*IJ1(N)^C3(N)*71(L) 
lji(.)*Ul(H  ♦:3C«)*71(N) 

02  (N)  s'J2  (N)  ♦A3  (  N)  *7?  (L) 

'J2  ( 1 1  »  02  ( L)  ♦  ft  3  (  N )  •  7  2  (  N ) 
CD>(3.«C.) 

PO  Iftfi  N=1*N® 

ii»r  00*00-01  (N)  *71  (N)  ♦J?  (  N)*V2(  N) 

rE7»5E/;o 
E7*E7^0E7 
CO*l./CSQT{T  (CD) 

x'ftos  (NEAL  (C)  )*A'5(AI»'AS(C)) 
:F(I<.fQ.l)  so, TO  3CC 
I*  (X,LT.E»r)  T<*1 
23 '  CONTINUE 
Pi.iNT  2120 

n  0» 


128 


20 


•  ?5 


70  r 


^5 


-*0 


<*•9 


50 


•*  ^ 


r,0 


71C 


•^le- 


37? 
71.  »• 


U4? 

45  0 

74  2 
40 


f  (<)=£/ 

FOLO(K)=EV 
X  =  A5r.(.^£AL<E\^)  ) 

Y=A3S ( AIMA3 (iV) ) 

Z~X*Y 

:F(Y/Z.wT.EF32)  L3(K)=1 
.3(0=2 

no  310  M=1,N? 

Vl(N)  =  i/l(M)*33 
V2(^()  =  y2(N)  ’OD 
>^(N,  <)  =t/l  (N) 
i;i(N)=  OKU)*-;? 

U2(*J)=  U2(N)*C3 

IFt-DCK)  .NE.l)  GO  TD  316 
x=(»(l(l)**2  <-1/1(2)  **2)*FV 
.F(X  .GT.O.)  GOTO  316 
t:(o=-E(<) 

DO  315  N=2.N»,2 
V(Ny<)=>y/(N,K) 

IF  (K.EQ.C)  GOTO  335 
CE»(0.,C.) 

CO  33:  N=l,N3 

C£*C£-^1  ( N)  *J1  (  N)  ♦''2  (M)  •ug  ( Ni 
TrJSo  >  ♦AOr  ( AI^AG  (OF) ) 

22:e,  K,  IT  ,£  (K)  ,  DEi/,  X 
I^<?.LT;|FS>  33  '3  735 
FF.I.JT  2u2l 
•'TOP 

00  340  N=1,N3 
ri(M)=^l(N)  ♦i/l('J) 

^2(M>=.<2((J)  <1/2  (N) 

.F(Xl.NE.O,C)  Sn  TO  45C 
J=L3(K) 

call  0^TH03(J) 

00  445  N=1,n5 
W1(M)=S1(N) 

V2(N»=S2(N) 

C  *  C  ♦  4 

IF (l5( K) . Nf . 3 »  MC=KO-? 

.■F(M3.EQ.nm?)  go  to  ICr 
C- 0  TO  4 : j 

c=o. 

:'»o. 

CO  342  N»2,N3,2 

CsC-Ul (H) *51 ( M) ♦J2(N)*F2(N) 
Ijb3«'J1  (S)  *?1(N>  ♦J?('<)  *r2  (M) 

C«5*2. 

0*0*2. 

DO  392  N»2«49,2 

-C’VKi) 

‘^l(N»  »S1(N)  -5*»/l  (N) 
'^2(M)»S2(M)-3*V2(1) 

‘'2(N)  =92  (N»  -:*rf2  (N) 

Vl(4»*3i(M) 

Yl('4)=3l(N) 
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V2(N)sS2(N) 

I’D  345  N=1,N3 
VI  (N)sSl (N) 

V2(N) =S2 (N) 

:"(L0(<)*I'1.Er3.T)  I^*sl 

IFCIC.Ea.N?)  GO  .'3  350 

GO  TO  430 

XSO. 

C— (3 

'JO  510  N=1«N» 
XsAI‘16G(E(N)) 

:F(X.L£.3  .;•)  33  T  3  60C 
£  I 'I)  =-£(N) 

C-“C 

no  610  J=2,N3,2 
V (J,N) sv ( J, N) *Z 
1  =  1 

IFd.GT.N?)  GO  TO  55C 
X=0. 

JO  540  .=I«N3 
fiC=-^.£:CL(E(l>) 

3F(X.EQ.3.>  «=AC 

IE  (A3. 31  .X>  33  ’’3  54Q 

J=i. 

Xsft" 

CONfiNUi 

c=E(:) 

Ed)  =E  ( J) 

E(  J)  =C 

JO  535  l=1.N» 

C=V(L,I) 

V(L.I)=V(L, J) 

V(L, J) =C 
I»I*l 
GO  TO  533 
C0'1’’IN'J£ 

<=3 

no  666  1=1, 

tr!=?EAL<Ed)) 
TF(AO.LT.O.C»  30  ro  666 
<=<♦1 
CQ'ITINUE 
I\*< 

ip»n 

!K»0 

00  560  Isl.NS 
a:  =  5EAL(£d)) 
]F(A3,LI.C.0)  go  to  62C 

‘■ld'4>=£d» 

00  536 

w(u,:n»*v(.,i) 

GO  TO  56o 

iF«!‘»*l 

•'l(I»)*Ed) 
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ro  557  L=1,N3 

55r  i>'(L,IP>=\/(L,n 

551  CCNTINIIE 

^0  555  1=1, N» 
rC555  J=1,NP 
tf(J,I)»W(J,l) 

555  P(X)*51(I) 

'  4^X,2^^  DU'PJT  •»)‘l  FF:£tfa  ) 

2  rXfjHllE?,,  5X,  PA?T, 14X,1 jHIMftS.  Pfi?T,l3X,  9M:ra  »!?• 

BX,14H0!frrjr'G:'4ALlTY  /)  •  .  < 

i-’ivp  sfi  ?ssi’4  ;;^6r?Mis;s;=}?:?5rsmi-?iE  ..i-i..  v--,,  .> 

-5T%‘?Li£feT'f3'S!." 

202"  FO^IIT  (43Hi  "AIL'JPE  TO  CONi/FiSF  2D  IT- RATIOS  ^T^P-  I 

TO  riMrj  F'SENi/E CTOS.  C^HOiO^Al  TO  THE  DT.'r 

1?"  ENO' 
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..  ,92(3:>  ,3!*  (33)  , 
a(3J) .  ^2(3:),  Ji  ( 


.  1 


15 


^3 


no  132 
- 

J=N+MO 
3  (N) 

20  55  .  =  JtM»ND 

c=a3(L)/Ae(M) 

i=i*i 

CO  3C  <=L»K.ND 

AB( Jl)sAB( Jl)  -0*A3(<) 

32  J1=J1*MD 

5?  vlcil  =V2(I) -:*'/2«'i) 
lO'"  \(2  (N)  »(/2  (N)  /  A3(  *J) 

».-S»|9 

V2(r))sV2(N)  /AB(N) 
no  250  K=1.N1 
tJ=N-l 
I=N 

‘•=‘175  ('^) 

CD  250  u  =  J.><.N3 

’  =  !♦! 

’5"  i2(N)=72(N) -A3('.)  *72(1) 

£Nn 


Mil 
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SU'i^'JUTIME  0<TH03(J> 

•  f  ^  M  ^  ^  ^  m 


:  U  \  ^  W  :  X  >  I  C.  U  N  I  n  W*  «>  \  # 

5S3l;f;  l!‘?5?r;hJ;l^:sy!!S[;'i?'-'-*^-^'-== 

COM-IDN  N0,?4P,N1  ,'l?,M3,CHS,Ik,  ir.TI,  AK  30)  «A3(3D),32(3e 
30  ,C3  (SO) ,  £11(3?)  ,'2(30  ,L0  (iQ  )  .>1rf3(3T  )  ,E  (3D)  .?1(3J)  ,  9 

#■*‘11  .  «9fT>l  .rni  nf^ii  .  u  it ;  _  i  fi  ti  n  t  r  \ 


ir 


3,J2(30) .  Vl(23) .  V2(3:)  ,FOLO(  3))  ,  »/(3i,,3C  )  ,A3(12C) 
C4=0. 

C2sD  • 

CO  ij  '(=2»NP,2 
f=M-l 

r,4=C4-Ui<M)  *31  ('•)  ♦J2  (N)  ♦]2(N) 

C2=:2-U1(U)  •?l(M)»J2(M)  *;2(^) 

"4=C4*2. 


32(30  » 34  (3  3)  «: 
92(33), Jl(3 


C2=02*2. 

:"(j.NE.o)  GD  ro  ZF 
no  23  N=2,NP, 2 
f*=N-l 

)'  =  ~4*V1(M) 
Fl(*1)=5l(M)-X-X 
X=32*Vl(N) 
n  (N) *31 (N) -x-x 
X  =  32*V2(‘1) 

‘'2(M)  =S2(M)  -X-X 
X*34*V2(N) 

2:  :2C4)*32(N) -X-X 
r,0  TD  4D 

2F  I'O  33  Ns2,NF,2 
V  =  N-1 

;  1('1)=31(W) -:4*Vl(M) 
M('l)*Sl(N)-32*\/l(M) 
'2(*()=32(«)-:2*V2(*() 
3'  :2(N)=32(N) -:4*v/2  (N) 
4r  ENO 


jiJj-.:CUTIi;i  V.r\JP  U,'i) 

IPL  <  r  (■)-..)  ,V  5  ) 

..iL 

^'=  N*r.' 

I  =  j 

L  =  » 

0  3  il 
■  =  1*1 

r'(  I):V(L,J) 
i-(L.LT.;<)  GO  Tj  12 
5-JtI 
F  0 i  T  (15) 
oJ  1 C  J=l, <S 
-  (JX 

w- ( JjsaiMiGCPC J)  ) 

:  i.jTIriuj 

OiwL  lM;-R7C(-<=‘,C",.4,=ri,OFI) 

00  2,  J  =  l, 

;  (  J)  =  C'lP!.X(-.FI(J)  ,CFT(J)  ) 

I  -  j 

i.  =  L  •*■! 

03  2,=.  J=1,  J 
1  =  I  ♦  1 

«/(u  ,  J)=F(H 

go  to  25 
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:  j j-ccu'^l.n  liivirT:  (AA:<,aai  .NP-*  1' I' 

AA^d  ,1)  ,  AmI  (1,1)  »Ar(1  ,  1)  ,hT(  1,  1) 

\.'Aw  A lA  K  ?,M!dXi  ,  <•<,  <I 
•TIUNSION  liPAC: 

'  '.'.<SNP*MP 

OJ  1=  1=1, N3 
!S=A^::(3,I)  =  D 
j.'J  !'■.  J=1,.JN,N» 

-■^lI,J>  =  AAx(I,  J) 

-:(I,J)iAAI(I, J> 

Ij 

iv  ^.TNTIriUi. 

■).)  lA^  1  =  1, U? 

-.'dXPsO. 

,••1X1  =  3. 
x^=C  . 

jJ  cw  J=1,Nk 

aP(  i:;pAc- (3,j) -1)  2n,oj,2: 
iJ  K<=.' 

■'j  S,  <=1,.JN,‘J  = 

<<=<<♦1 

oO  -3  <!=<<,<< 

IP(I-rAC>(3,<l)-l) 

3:  -.or^TiMOv 

j  ,■!  : ■( mP (j,<) ♦Mn.t j,K) ♦a: ( j,^) -uT ( J, K) -X?)  i3,‘-‘;,io 

•»>  irHWsJ 

i:c=K< 

-  -IAXksA^C  J  ,<) 

1  i-XisAI ( J,K) 

<'  =  A<-.AX\*A  ^A  X-<-H;-U  <I*A!!A  Xi 
•i3  :w.nTI‘!U£ 
oa  CjlHlN'J^ 

I  J0LU**sICC 

Ii3*,Ci:(3,I«0LU'1)  =iaPAw^C*,ICOL'.r'!)»  1 
>  IP ( ir.j.^-icoL'j'-i)  r;,'s3,7: 

'O  .■<.  L=l,iJN,-4P 
X  ^=Ar.(:.<OW,^) 

''Is  mT  (I-XOW.L) 

-  '.( I  <OW,L)  =AF  (X:aLU~!,L) 

.  II ( 'OK,i.)  =Hi(:cot  JM,L) 

■*i  llCOuJM,-)  =  x» 

5.’  i:  ( ICOLJ'<,L)  =  XI 

a:  .jFlat  (1,1)  sl^iOW 

l5  =  AL‘^(2,I)  =  IC0aU'1 

I  ;=  ( icouu’''-i)  ••(■■>♦1 

!•<(  laJlUMtlw)  =1.  J 

-.iriC0LJ'1,Iu)  =  G.  3 

30  lal  L=l,fif4,*JP 

X  •>=  A.lAX.t*A  -!AX  =  +  AMaXI*  Ar-'AXI 

a  XlsAF 

A‘.  (ICOLJMfA)  =  (Xl*A>lAxp4AI(i;aLUM,L)  •A''Avl)  /x  = 

AI  (  IC0LU.-I,L)=  («I(IC0LJ.i,i.)*AflAX-;-Xl*f''AXi)/XO 

131 

Zj  IZz.  L1=1,NS 
Tr(Ll-I3JLU1)  lli;,l’3,a.l. 

11a  X  5=  An.(Ll  ,i:» 

Xil-l(>l,Ia) 
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*  w 

13C 

1 


—  I —  J«C 

-  ;(Llil>>  =  -•.* 

jT  L  = 

CLXtD  — Axfcl 
:jN*1:IO:; 

'  jljllN  J" 
Zi'*'’lUUZ 


l-i'CirOLOM,  l1  C1C''..J?< 


17.  1=1, !C? 

L  =  '1P»1-I 

I*(iS-A-iCl,i.»-TS’43=C?,L»)  i5.,37-,l»- 
j  <0  J<=  IS^-  3i  C 1 ,  L» 

J3:s!.U-=xS^-iC£  C2,.» 

•»;?♦! 

J3=  IJ~J!.U?'.-l>*N3*t 
OJ  IrS  <=l,aP 
ir-:=i-iCK,  J?l 
<r=Aic<,  j-.i 

arcii,j?»=Aii<,  j.) 

1t.C<,JC»  =  <i? 

ZJIZVSJ- 


L)-*T 
•.»  ’Jf-C 


9 
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1 


1? 


-a 


^3 


.1 


5' 


;U=»’DUTiNE  MDDft. 

CO'iP-tX  P,Et-n,9l,R2.l«JjU2j(/  I.WEW.A3,C 

Ai  j72» « a3f  35)  «5?J^£l ? 

2?t).  23  (22>.i:i(33)i'2(3n.Lr(>Q)«'W  3  (2  j).v<  30. ^113?).  P2  132),  Jl( 
3,U2(32)  .  »/l  (33)  ,  V2  (30  ,  FOLD  (  31)  ,V(3i, 30, 03(120 

Kr-<c 


?:  N=M*i 

C=(C.,-1.) 

I=»3(N) ♦! 

Ir (1-2)13,21, Au 
ir  5i(=tI»i4G(E(K) ) 

l'^(Erf.Ll  .C.)  GOTO  12 
r (S)=C0\J6(E(V) ) 

'»C  11  J=1,NF 
11  V(J,'i)*20NJG(V(  J,>l)» 
1?*  E(K)  S-2GNJ3  (£  (X) ) 
LC(X»=3 

13  j=i,MP 

1?  V(J,<):COriJG(tf  (  J.X)) 
to  1%  J=2,SC,2 
10  V(J,<)=-tf(J.<)»r 


<=<-1 
GOTO  53 

?«■  ?V=^££;.  (£(•(»» 

£(M)  sEV 

>f=4?GI^£4L(V(l,f<)  ))  *A?r(=-cAl  (  /  (2.>ll  )) 
lP|y.5T,l.£-l2l  G0~0  3*- 
»C  33  l.=l,'IP 
3’  iflL.')l=*I"AoH»(;..'l)» 
bC-3  53 

?fc  ro  35  L=1.»IF 
2f 

ro’o  52 

v:  £y=«l*(AG(£(N)) 

OF  (£¥.1.7.3.)  Sr  0  fcl 
£¥=-£¥ 


o:  E|:«=2'(».*(t.,£¥l 
?:  ro  51  J=2.¥P.2 
51  ¥IJ,*ils¥IJ,M  *3 

:^(.ocii .K£.2i  srro  i?*^ 

X=A3>l^iLt¥tl«¥)»)*A5ri^£lL;  ¥(£•¥))) 
IFI«.eT.l.£-l?)  GOTO  5*. 

"♦Cf  53  Csl.HP 
57  ¥(t.,¥lsAlNA&(¥(L.¥)* 

GOTO  IJi 
Zi  55  :.sl.!¥F 
St  ¥  (L*^)  «%£«!(<#(!.•  91)  > 

131  ri*G 
•*-»r 

llT 

lF(:.n  j)  .to.') 

IF(_3(  j).£a,?) 

17-  C0^1’»^J£ 


li 


I  »r»> 
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o3 


■■.5 


7? 


71; 


It) 


'  5 


::=N»-ic-ir 

1*1 

ll'*  IFCI.GT.IK)  SPTr  133 
X  =  3. 

np  120  LsI.N? 

lF(.D(i.)  .Ni.l)  COrO  120 

CC*A9r;(xEa-(£(L))) 

IFtX.EQ.i.J  <*4C 

IF  (flC.GT.X)  GOTO  12D 

J='_ 

X  =  fl2 

12''  CONTINMi 
l=-3(I) 

LO(I)sLD( J» 

■.0(J»  *L 
r=£(I) 

E  C)  *£  ( J) 

E(  J)  =C 

LO  125  .=1,N» 

C=\/(L»I) 

»/(i.,r»=v(L.  J) 


1?? 

13'“ 

y(,. j)=: 

:=iii 

GCfC  lie 

(I  .C-T.Ir  +  IC) 

core 

15C 

’5 

x  =  3. 

1 0  lUO  .=I,N3 
iF(L J(L) .NE.3) 

GOTO 

14P 

aC*’£AL(E(L)) 
!-(x.EQ.u.i  <=a: 
iF(a:.;f.x)  goto  luo 
j*» 

x=a: 

1!*:  continue 

L=LD(T) 

Lr(I)=LD(J) 

LO( J) =. 
rsE(!) 

F (!) *£( J) 

E(  J)  =C 

•'O  135  l  =  1.N» 
C=y(L,I) 
t/(LtI)*i/(L.  J) 

135  V(LfJ)*C 
]=I*1 
r.rTp  133 

157  :f(f.g£,NF)  GOT')  17" 
Xsl, 

PO  150 

AC*A3f {AIMAG<E(L)  )  ) 
lf(x.EQ.o.)  xssr 
IF(A-,g1  .X)  G'lTn  16" 
J*'. 

Xs  AC 

15"  OONTINUE 
L»l.O(I) 

I  0(1)  *i.3(  J) 
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!  t 


:t3 


’5 


-0 


U5 


=50 


-5 


•■in^TUTIIIE  JTIsr  >(y,><1,G1,»H,  AO 

i-I’iEMSION  V  (4«?)  «  E(3>  «AJ(2,  £1  ,  Z  (?t  4)«  3  (4)  «Vj(3)  »  A(3t 
1,A  K4,4J  .ArHS.*. ) 

COis.EX  Q(3.3).CA{3,8),AK(3,4) ,Ar(5,3, 9),P(3) 
T(l)=-.77459=)53n24148 
t (2) =C. 

W(l) =.55555555555355 
W(2) =. 355666899 
£(3) =-E(l) 

W(3):M(1) 

L  =  1 

CO  4  ZA=1,3 
no  4  J As  1,3 
Z(1,1)»-1^£(IA) 

?(l,2)=l-E(Ifl) 

2(1.3)  =1«-£(1A) 

2{l,-*)=-l-E  (I  A) 

Z(2,l)=-l*E(JA) 

Z(2,2)s-1-£(JA) 

2(2. J)=l»£( JA) 

2(2,4)sl-£( JA) 

:  (1) sZ (2.1) *2 (1 .1)74. 

=(2)=-Z(2.2)*2(1.2)74. 

?(3) -Z (2.3) *2 (1,3)74. 

?(4)=-Z(2,4)*Z(1.4)/4. 

•10  1  1=1.2 
:o  1  j=i.2 

AJC.  J)SJ. 

50  1  K=1.4 

1  AJd  ,  J)=AJ(I,  J)  ♦2(I,<)*Y(K.  Jt  7  4. 

0E'!'sAJ(1.1)  •AJ(2.2)-AJ(1.2)*1J(2.1) 
fJl=AJ(l.l) 

AJ(1,1)=AJ(2,2) 7jET 

AJ(1.2)»-AJ  «1,2)75eT 

AJ(2.1)=-AJ(2,l)/3Er 

A J(2.2)=AJ170ET 

'lO  2  J»1.4 

JJ=J4J 

.0  2  3  I  si. 2 
A(I, JJ-1)=3. 

?*  A(I.JJ)»0. 

**0  2X  K*X  ? 

A(l,  JJ-I)2a(1,  JJ-1)4AJ(1,K)  •!(K,  J)7  4. 

’1  A(2. JJ)=4(2,JJ) ♦AJ(2,K)*Z(K,J) /4. 

A(3. JJ-1)=M(2, JJ) 

?  A (3.  JJ) sA (1, JJ-1) 

DC  11  1=1,3 
"0  11  J«l,3 
11  C(i,j)=b. 

C(l.l)»Hlt5l*'l 
C(2.  2)*rtl^3lfr,l 
C(1,2)*M1 
C(2,l) *rtl 
C(3, 3)»G1 
'0  2  1=1.3 
no  3  Jsi.fc 
CAd,  J)s!3. 
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-1 


'5 


*0 


3r 


a 

ir 

I. 


''O  3  Ksi,3 

CA(I  ,J)sCA(I,  J)  «^C(I«<)  •&(<«  Ji 
"D  ?]  I=lt3 
OC  33  J=I,5 
'K(r#j)=3. 

Lrf C  3  3  K®  *  3 

Hk(I.  J)sSk(1,  j)  «A  (K«!  )«C«  «•  I) 

•:o  e  1*1. <♦ 

rj  9  j»it4 

Aid,  J)=WF*c*(*3  (U*9(  J) 

'0  9  1*1.4 
OC  9  J»i,4 
:i*i »i 

JJrJ»J 

AN(i:-l, JJ-1) SA1(», J» 

AN(II» JJ) J) 

rc  13  1*1, b 

ro  11  Js!,9 

t  »!.*1 

r.c  5  1*1, e 

5  j=i,5 
nc.  7  <*1,3 
» (<)  =  ?. 

■_*3*  (<-l) 

no  7  M»l,3 


t  (<)  *P  (<>♦*<  (1)  •fir  (j,  j,  ^,Q) 

AKd.j)*:. 

CO  5  1=1,3 

A<d,J)sAK(!,J»4H(M)*Pf1) 

FND 
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5:'J3’3Uri:'lE  W^D5C(LC0JK‘<  .NS;) 
Cr'11DN/3ANft?.3/M‘I,'ir.NL)*’^LK.  S(o»i )  .ft  (o9. 
CU'i».EX  4,9.5 
r4U'13;,K*N'JM3L<»l 

(2)  (Jill)  ,  {4(11,  Jl)  ,Jlsl,—.)  ,11* 
3i;  11=1, NN 
LC3UNT*i.C0l)NT»l 

if(l:ount-ns) 32r.32: ,3?o 

3?'"  *,Z*N*I*I1 

EI!1)*3INZ) 

E IN?)s2. 

r.j  313  jisi,^^ 

{ (II. J1)>A(NZ, Jl) 

•'1'  '“•(NZ,  Jl)='.  . 

■«3r  E'.n 


34)  ,r(3C.3:) 
1  ,NN) 
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